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FINITE ELEMENT ANALYSIS 


OF LOW SPEED VISCOUS AND INVISCID 
AERODYNAMIC FLOWS 
By 

A.J. Baker and P.D. Manhardt 
Computational Mechanics Consultants 
Knoxville, TN 37920 


SUMMARY 


A weak-interaction solution algorithm is established for predicting 
aerodynamic flow about isolated elementary airfoils. The procedure employs 
finite element numerical methodology applied to solution of the partial 
differential equations governing isentropic potential flow, and the viscous 
and turbulent boundary layer and wake flows downstream of the sharp trailing 
edge, Iteration within the algorithm accounts for computed viscous displace- 
ment effects on the potential flow, the modifications to which in turn alter 
the viscous flow computations through imposed freestream pressure gradients. 
Closure for turbulence phenomena is accomplished using both first and second 
order models. 

The C0M0C finite element fluid mechanics computer program was modified 
to sequentially solve the identified equation systems for two-dimensional 
flows. A comprehensive numerical experimentation program was completed to 
determine factors affecting solution accuracy, convergence and stability for 
the combined potential, boundary layer, and parabolic Navier-Stokes equation 
systems. Good accuracy and convergence are demonstrated for potential flow 
prediction, including angle of attack and a curved wake trajectory. An eco- 
nomically accurate discretization procedure is identified for turbulent 
boundary layer computations which provides adequate sub-layer definition 
coupled with good solution speed. The parabolic Navier-Stokes equations, 
with turbulence kinetic energy closure are shown to accurately predict near- 
wake turbulent mixing processes for dissimilar turbulent boundary layer merg- 
ing downstream of the trailing edge. Each of these solutions is obtained 
within the identical finite element framework of C0M0C. In combination, with 
iteration through a sequential solution, the algorithm can establish an accu- 
rate and stable procedure for analytical characterization of complete flow- 
fields about practical elementary airfoil configurations. The results of 
this evaluation, which document existence of a potentially viable supplement 
to wind tunnel testing, are presented in this report. 


INTRODUCTION 


Aerodynamics is a unique branch within engineering fluid mechanics. 
Elsewhere, the non-linearity introduced by the convective acceleration 
term in the Euler or Navier-Stokes equations precludes closed-form solutions 
in all but the most elementary cases. In aerodynamics however, for all low 
Mach number flight and slender bodies up to transonic Mach number, the differ- 
ential equation governing isentropic flows is the familiar linear elliptic 
Laplacian. As a result, potential flow aerodynamics developed rapidly and 
became a classic study in engineering mechanics (cf., ref. 1). Concurrently , 
it was equally well recognized that potential flow analysis is only a first 
approximation, with the major limitation being primarily neglect of viscous 
and wake effects. Furthermore, accurate determinations for high performance 
lift systems, including leading edge slats and trailing edge slotted flaps, 
is well beyond the elementary analysis, with the result that the wind tunnel 
has been universally employed to accomplish design optimization. 


The advent and ready accessibility of the large scale digital computer 
in the mid-1960's, and the phenomenal growth in its operating capabilities 
since then, initiated an entirely new approach to aerodynamics analysis. 

While the "computational wind tunnel" may still be a distant vision (cf., 
ref. 2), the use of numerical solution of the governing differential equations 
has emerged as a viable alternative to exhaustive experimentation in some 
cases. For example, a broad class in potential flow is now amenable to 
solution using numerical techniques to evaluate the singularity distributions 
imbedded within the Green's function analytical solution of the governing 
Laplacian equation. Specific elementary flow singularity functions, e.g., 
source, sink, vortex, doublet, etc., are combined to mathematically model 
thickness, lift, and the flow tangency boundary condition. Hess and Smith 
were among the pioneers in the development of the concept and its early re- 
duction to practice (e.g., ref. 3). Refinements have extended applicability 
to supersonic flight as well. The proceedings of two recent conferences, 
held at NASA Langley Research Center (ref. 4,5), document the tremendous growth 
in this area and the versatility and utility of the concept. 

However, as mentioned, potential flow determination is only a first 
approximation to the solution and must be corrected for viscous effects. The 
concept of the weak-interaction solution algorithm assumes that viscous (and 
turbulent) flow effects can be included by determination and use of an "effec- 
tive" inviscid shape. These corrections are established by determination of 
the boundary layer character, and everywhere augmenting the airfoil thickness 
by the computed distribution of boundary layer displacement thickness. The 
elementary integral form of the boundary layer equation is typically utilized 
(ref. 6), with an experimental correlation function for shape factor. This 
concept has proven quite useful for elementary two-dimensional airfoils at 
modest angles of attack, and has been extended to analysis of high lift sub- 
sonic airfoil systems including an accounting for merged slot flows (ref. 7). 
However, since the viscous flow methods are elementary, at best, there exists 
ample area for improvement of the accuracy of the analysis. In particular, 
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as noted by Morgan (ref. 4, p. 730), the state-of-the-art, viscous-augmented 
potential flow solutions are amenable to improvement in determination of 
profile drag from downstream wake characteristics, flow characteristics after 
boundary layer separation, and pressure corrections in slot areas between over- 
lapping elements of multi -element airfoils. In each case, the improvement 
required is in the handling of the viscous flows associated with practical 
airfoil configurations. Equally important is improving the stability of the 
potential flow algorithm when viscous corrections are included in an iterative 
analysis. 

The mission of the present study was to examine the viability of the 
finite element numerical solution technique applied to the differential equa- 
tion systems governing weak-interaction theory in low speed aerodynamics. 

Initial impetus was provided by the observation that the finite element algo- 
rithm is equally applicable to the potential flow and the turbulent boundary 
layer and wake flow governing differential equations. Furthermore, since 
finite element methodology is very flexible regarding discretization non- 
uniformity and gradient boundary condition specification on arbitrarily orien- 
ted (curved) surfaces, it appeared intrinsicly well suited to the aerodynamics 
problem class. In particular, since the flow tangency boundary condition can 
be applied directly on the "effective" inviscid contour, including the wake 
contour predicted by a complete viscous and turbulent flow solution, it 
appeared that iterative stability should be acceptable within the limits of 
validity of weak-interaction theory. 

The results of this study, as generated using the C0M0C finite element 
computer program, are reported herein. Although ample results on numerical 
experimentation have been reported on finite element solution for linear 
potential flows (ref. 8,9), including airfoil configurations (ref. 10-15), a 
detailed analysis of the aerodynamic accuracy has not been documented. Factors 
affecting solution accuracy, including use of highly non-uniform computational 
meshes, and reduction of finite element results to equivalent pressure dis- 
tributions, have been thoroughly evaluated. The capability to automatically 
generate required discretizations and efficiently cycle through angle of 
attack are now documented. Finite element solution of the attached viscous 
and turbulent flows and the wake flow downstream of a sharp trailing edge 
of an airfoil have not been previously reported. This combination constitutes 
the weak-interaction problem and is sketched in Figure 1. The complete viscous 
flow equations governing these additional regimes in aerodynamics are developed, 
including models for both elementary and higher-order closure for turbulence 
phenomena. Numerical solutions as generated by C0M0C, document an acceptable 
accuracy attainable with a minimum fineness discretization that is highly non- 
uniform for both boundary layer and elementary wake flows. The combination 
of these solution capabilities, existent within a single general purpose com- 
puter program, is then indicated to provide a convergent weak-interaction 
solution algorithm for a representative isolated two-dimensional airfoil 
configuration. 
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SYMBOLS 

boundary condition coefficient 

Van Driest damping function; constant 

body force 

isentropic sound speed 
coefficient; chord 
drag coefficient 
skin friction coefficient 
lift coefficient 
differential 

finite element solution energy 
function of known argument 
Froude Number 
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metric coefficient; mesh parameter 

boundary layer shape factor 

thermal conductivity; turbulence kinetic energy 

generalized diffusion coefficient 

differential operator; length scale 

differential operator; length 

finite element index 

Mach number; number of finite elements spanning R 
unit normal vector; nodes per element; dimensional ity 
pressure; generalized parameter; convergence parameter 
Prandtl number 

generalized dependent variable 
generalized discretized dependent variable 

universal gas constant; domain of elliptic differential operator 

Reynolds number 

coordinate tangent to curve 

finite element assembly operator 

time 

temperature 
velocity vector 
reference velocity 
friction velocity 
scale velocity 
molecular weight 
Cartesian coordinate system 



friction velocity Reynolds number 
angle of attack; parameter 
ratio of specific heats 
circulation 

closure of solution domain R 
Kronecker delta; boundary layer thickness 
boundary layer displacement thickness 
increment 

turbulence dissipation function; thickness 

boundary layer momentum thickness 

transformed coordinate 

transformed coordinate 

Karman coefficient (MLT) 

multiplier; turbulence sublayer constant; 

dynamic viscosity 

kinematic viscosity 

density 

mean flow Stokes stress tensor 

Reynolds stress tensor; wall shear; 
integration kernel 

perturbation potential function; finite element functional 

total potential function; 

generalized initial -value coordinate 

turbulence damping factor 

global solution domain 



Superscripts : 


e effective value 

T matrix transpose 

+ turbulence correlation function 

mass-weighted time-average 

time averaged 

A unit vector 

' mass-weighted fluctuating component; ordinary derivative 

* approximation 


Subscripts: 

°° global reference condition 

e freestream reference condition 

i,j,k,4 tensor indices 

non-tensor index 

m finite element domain 

n normal 

o initial condition; stagnation reference 

t time derivative; turbulent; tangential 

w wall reference condition 


Notation: 


{ > 
[] 
u 
n 


col umn 
square 


matrix 

matrix 


£ 


union 

intersection 
belongs to 
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THE FIELD EQUATIONS IN AERODYNAMICS 


The description of a state point in multi -dimensional aerodynamics is 
contained within the solution of a coupled nonlinear partial differential 
equation system describing conservation of mass, momentum, and energy. Unique 
solutions are obtained upon closure of the system by specification of consti- 
tutive behavior and boundary conditions. In Cartesian tensor notation, the 
conservation form of the governing Navier-Stokes equation system is 


L(p) = If + afr<P u j) = 0 

J 

Lfpu j) = & pu j> + 3xr[ pu i u j + p6 ij 

0 



+ pb. = 0 


L(pH) 



a. -u. 
U i 




(1) 

( 2 ) 

(3) 


The dependent variables in equations ( 1 ) - ( 3 ) have their usual interpre- 
tation in fluid mechanics where p is mass density, Uj is the velocity vector, 

p is the static pressure, b is the applicable body force, and H is the 
stagnation enthalpy. For a single species fluid, 


H e 



dT + 2 1 Vj 


(4) 


where c p is specific heat and T is static temperature. Furthermore, q. is the 

heat flux vector and a., is the Stokes stress tensor, defined as 

13 







2y 3u k 
3 9x k °ij 


(5) 

( 6 ) 


In eouations (5)-(6), k is the thermal conductivity and y is the dynamic 
viscosity. An equation of state closes the definition; for example, for a 
perfect gas 


p » pRW _1 T 


(7) 
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where R is the universal gas constant and W is the molecular weight. 
The important non-dimensional parameters in aerodynamics are: 


Reynolds Number: 

_ ^00 00 
Re = 

^oo 

(8) 

Prandtl Number: 

D - C P li 
Pr =k 

(9) 

Mach Number: 

8 3 

IK 

«ls c 

(10) 

Froude Number: 

P _ u l 
Fr = gl 

(11) 


In equations (8) and (11), l is a scale length of the problem and c is the 
isentropic sound speed. In non-dimensional form, the Navier-Stokes system 
becomes 


+ pS-jj - Re" 1 a. J + Fr" 1 pb i = 0 

(y- a -jj u i - Re ' 1 Pr _1 p ~- 

. (14) 

- (y-i)M^g|= o 

Equations (12)- (14) form the basic system governing all classes of aero- 
dynamic flows. Their straightforward solution represents a formidable task, 
and is not required for many flows of practical interest. The primary reauire- 
ment is to determine the pressure distribution produced by flow about an aero- 
dymanic shape, and to evaluate the resultant induced pressure and viscous 
drag forces. The concept of weak-interaction assumes existence of a conver- 
gent iteration procedure between a potential flow determination about the ' 
effective inviscid airfoil shape, and a viscous and turbulent flow solution 
in the immediate vicinity of the airfoil and in the wake. The weak-inter- 
action hypothesis is usefully valid at modest angles of attack, wherein the 
viscous flow is parabolic, i.e., the governing equations obey a generalized 
boundary layer order-of-magnitude simplification. In this instance, the 
inviscid solution pressure distribution can be applied unaltered to the viscous 
flow solution. In turn, the displacement surface shape of the latter determines 


L(p) E at + fcnW = 0 

J 

L(pu i ) s^tpu,) + 

L(PH) E £(PH) ♦ affpUjH - 


(12) 

(13) 
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the effective inviscid airfoil contour, and thus modifies the computed in- 
viscid pressure distribution. For these flows, the iteration between in- 
viscid and viscous flows should converge to a unique solution. The remainder 
of this section develops the various required governing differential equation 
systems including closure for turbulence. 


Aerodynamic Potential Flow 

The Navier-Stokes eauations (12)-(14) are simplified to describe the 
steady, isentopic irrotational flow of an inviscid fluid. The Stokes stress 
tensor a-jj vanishes identically, and H equal to a constant is the solution 
to equation (14). The potential flow description is then 


L(p) 



pu j) = 0 


L(pUj) 


3x j- 


p Vj 


P6 t 


= 0 


(15) 

(16) 


These (four) equations can be usefully combined into a single equation using 
the definition of isentropic sound speed 


c 


2 = 


V 

dpjs 


(17) 


Expanding equation (15) and replacing the resultant density derivative with 
equation (17) yields 


L( Uj ) 

For isoenergetic flow, where 
state 


n 9u_. 


- c ' 2 u i u j 


ax. 


j _ 


= o 


(18) 


subscript zero indicates the stagnation reference 


c 2 = c 2 - XlJ. u . u . 

c 0 2 u.u. 


(19) 


Equation (18) is a highly non-linear, first-order partial differential 
equation valid for all Mach numbers. As a function of M, it can display 
elliptic, parabolic and/or hyperbolic differential character. The class of 
interest corresponds to an irrotational (or rotation preserving) velocity 
fieldv ' Since the curl thus vanishes to within an arbitrary constant 
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( 20 ) 


u. = -~~ 
3 3x j 


where $ is termed the velocity (total) potential function 
tution, equation (18) becomes 


L(*) = 




- c 


“2 


3 $ 3 $ 


3x n - 3x 


j] 


3 2 $ 


3x i 3x j 


= o 


By direct substf 

(21) 


Equation (21), a second-order partial differential equation, preserves 
the mixed differential character cited previously. For subsonic flow at small 
Mach number, equation (21) simplifies to the elementary linear Laplacian 


iw * 


3 2 <j> 

3x..3x. 


= 0 


( 22 ) 


Boundary condition specification for equations (21) and (22) requires 
enforcement of flow tangency to given coordinate surfaces. Assuming nj the 
unit outward pointing surface normal vector, this boundary constraint becomes 
(from equation (20)) 


^ _ii 3$ ^ 

u j n j = U n 3x7 n j 


(23) 


Note that U n vanishing identically enforces flow tangency. In the instance 
of small Mach number equations (22)-(23) represent the classic von Neumann 
specification for an elliptic differential equation. A unique solution is 
obtained by setting $ equal to an arbitrary constant at one point in space. 


While total potential has been evaluated in finite element literature, 
the alternative perturbation function definition is common to aerodynamics 
and has been determined preferable for the finite element approach as well. 
Herein, the velocity field is referenced to the undisturbed freestream as 




(24) 


where £j is the freestream flow unit vector related in an elementary manner 
to angle of attack. Substituting equation (24) into (18) and dividing by 
the freestream sound speed 


L(*) = 


6,, M z e 
ij 00 




+ 


e. 

+ £ 

_ 1 

n 3x j 

j 3X 1 - 

3x- SXjJ_ 


3 2 4> _ 


(25) 


3x* 3x . 


= 0 
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Equation (19) has been referenced to freestream. 



+ 



3(fr 3(j>~ ] 
3x • 3x.J 

J J 


( 26 ) 


Sufficiently remote from the airfoil, since Uj = U^j , the appropriate 

freestream boundary condition is <j> a constant. However, the Kutta condition 
for angle of attack modifies the interpretation to allow a jump in 4> across 
the wake trajectory. Hence, the farfield constraint valid for all cases is 


$7"i = 0 < 27 > 

where n, is the outward-pointing unit normal ^vector. At the effective in- 
viscid-flow surface of the airfoil, since ujnj vanishes identically, the 
boundary condition on <p is 


3£ 



A 



e -n . 
J J 


(28) 


Hence, angle of attack and local airfoil contour (unit normal) are applied 
directly as a boundary condition. 


Equations (25)-(28) are valid for all Mach number and display the mixed 
differential behavior. The slender body approximation can significantly 
simplify equation (25) by the neglect of all products in perturbation velocity 
and the assumption that ej is aligned with a principal coordinate, say x 2 
i . e . , 


e . = 6 M 
J Jl 


Equation (25) then takes the form 

L(*) - - ^U«jl + (Y-l)M^ 5ij 


+ M“ 


' s njt + 


d Z <$> ^ 


9x i 3Xj 


0 


J v 

In expanded form for two-dimensional flow, equation (30) is 


(29) 


(30) 


L(4>) = 


1 - + (Y+1)M^ 



+ 2M 


2 3(f) 
°°3X 2 


3 2 <j) 


+ 


+ 


(y-i)m^ 


1 J 



(31) 
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which is recognized as the complete small disturbance equation (cf . , ref. 16). 
For most flow velocities below transonic, equation (30) or (31) can be line- 
arized by retention of only (1-M<|) on the lead term. The Prandtl -Glauert 
(ref. 1 ) coordinate transformation can then be employed to eliminate ex- 
plicit appearance of M^. Hence, for all subsonic flows about slender bodies 
equation (30) becomes the Laplacian, i.e., 

L <+>*33$&7 =0 ( 32 ) 

As mentioned, the desired output from the potential flow solution is 
the pressure distribution over the effective airfoil surface and onto the 
viscous trailing edge wake. The pressure coefficient is 



pu-u. 
H i i 


P U z 

oo oo 


(33) 


Substituting equation (24) yields 


r — g 3<j> A 9<j> 9<j) 
p 9x- e j " 3x. 3x . 

J J J 


(34) 


Evaluation of equation (34) is facilitated by resolution of the vector 
field into scalar components locally garallel and perpendicular to the 
effective airfoil surface. Letting £j denote the unit tangent vector, and 
(n,s) the corresponding surface-oriented coordinate system. 


= M-n + M-t 

ax, an j as c j 

J 


■ + 


(35) 


The last form was obtained using equation (28), Substituting equation (35) 
into (34) yields 


C = (e-n.) 2 + 2-^-e.t. - 
P j J. 3s e j\j as as 


(36) 


Recall that ej is the constant unit vector for U^; the remaining terms in 
equation (26) are all a function of the surface coordinate s. 
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The decomposition in equation (35) has yielded C D in a form directly dependent 
upon airfoil geometry. Hence, ft. and u are M potentially available to any 
required degree of accuracy. As will evolve, the finite element solution for 
distribution of <£(s) is a direct computational result for the zero normal wash 
boundary condition equation (28). Only these surface values are required to 
evaluate equation (36), using for example, spline interpolation techniques to 
establish the required tangent slope distribution. 


Aerodynamic Viscous Flows 

The complete Navier-Stokes system equations (12)-(14), is required to 
adequately describe general viscous flowfields including regions of recircula- 
tion. For small angle of attack, such that the viscous flow remains attached 
to the airfoil surface, the considerable simplication of boundary layer theory 
(e.g., ref. 6 ) is appropriate. These equations in general are an insuffi- 
cient description for the viscous wake flow downstream of the trailing edge 
terminus however, wherein the parabolic Navier-Stokes equations are required 
( at least ). The boundary layer equations are a sub-set of the parabolic 
system, the establishment of which is readily accomplished. 

The three-dimensional “parabolic Navier-Stokes" (3DPNS) equations des- 
cribe the steady time-averaged confined or unbounded viscous and turbulent 
flow of a compressible, nonisoenergetic fluid wherein: 

(1) A predominant mean flow direction is uniformly discernible, 

(2) In this direction (only) diffusion processes are negligible 
compared to convection, and, 

(3) No significant flowfield disturbances are propagated upstream 
against the predominant flow direction. 

To establish the system, assume the Reynolds decomposition of the instantaneous 
velocity vector^ u q * , equation (2), into a steady mean flow u-j and velocity 
fluctuations uf as (cf,, ref. 17) 


u-j = u. + ul (37) 

In equation (37), u-j is interpreted as the mass-weighted, time-averaged 
steady flow, i .e. , 


u 


i 


PU. 


p 


(38) 
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By definition, using the overbar to denote time-averaging 


t +T 


pu 


1 r° 

i - Hm j (pu 1 - pu. )dt = 0 
‘ I-*-'® i 1 


( 39 ) 


The concept of time-averaging yields the important relation 


pu 


• U • - pu .u • + DU -U • 
1 J J ^ 1 J 


(40) 


For the isoenergetic steady flows of present interest, the time-averaged 
Navier-Stokes equations for the mean flow become 


9(pu.) 

L(p) = -k=±-~ 0 


3x i 


KP a i ) P U j 3Xj + 3x^| ' 3x.hj ' pU i U j] " 0 


(41) 


(42) 


In equation (42), a-j -j is the time-averaged mean flow Stokes stress tensor, 
equation (6). The second term in the divergence is called the Reynolds str 
tensor, x-jj 


stress 


T • • = -PU • U • 
1J H 1 J 


(43) 


The parabolic approximation to equations (41 ) - (42) constrains appro- 
priate summation limits. Assuming the X! coordinate aligned with the direc- 
tion of predominant flow, this approximation becomes 



= _H_(i . 
Re u 






3x. 



(44) 


The mean flow unidirectionality also affects retained scalar components of 
the Reynolds stress, to be developed. The subscript bar notation denotes 
the index is not eligible for the summation convention, but merely takes on 
the value of the synonynous tensor index. In expanded form the 3DPNS system 
is 
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Note that equation (46) represents three scalar components. The two-dimen- 
sional system, i.e., 2DPNS, is obtained by deletion of all x 3 partial de- 
rivatives. The dependent variables of the 2D and 3DPNS equations are the 
steady-flow mean velocity vector q. and pressure p . The presented 
system is not closed since the scalar components of the Reynolds stress tensor 
remain to be identified. This is accomplished in the following section. 


Over the modest range of angle of attack in subsonic flow to which the 
present analysis is restricted, both the boundary layer and wake flow diff- 
erential equation systems are subsets of the 3DPNS system. For^the former, 
the first and third equations in (46) are solved for u 2 and u 3 respectively, 
as initial -boundary value statements. The mean velocity component normal to 

the surface, u 2 is determined from ini tial -val ue solution of equation (45) 
in terms of its prescribed value at the airfoil surface. Through an order 
of magnitude analysis (cf., Tennekes and Lumley, ref. 18), the u 2 momentum 
equation can be simplified to 


L(pu 2 ) 



pu 2 u 2 


= 0 


(47) 


Therefore, pressure variation through the boundary layer remains a second- 
order phenomenon, and is linearly proportional to a scalar component of the 
turbulence kinetic energy. Hence, to first order, the inviscid pressure dis- 
tribution determined from equation (36) is applied directly to solution of the 
viscous boundary layer flow. Equation (47) could be evaluated from frees tream 
where the pressure is known to the airfoil surface to provide a second order 
correction. This same order of magnitude analysis identifies -pu v ti£, 2<-^3 
as the dominant components of the Reynolds stress tensor. 1 * 5 

Determination of the viscous and turbulent flow in the wake downstream 
of the trailing edge using the boundary layer equations is inappropriate 
except for zero angle of attack and a symmetric airfoil with a cusped trailing 
edge. Solutions in the immediate vicinity of a blunt trailing edge would 
require the full Navier-Stokes equations. For a sharp trailing edge whereat 
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non-symmetric attached viscous-flows merge, the 2D or 3DPNS system can provide 
prediction of the ensuing wake. Guidance for the solution procedure comes 
from available boundary conditions; initial conditions for u-j and p are pro- 
vided by the individual boundary layer solutions at the trailing edge. Down- 
stream, Ui, u 3 and p are known as a function of Xi at both freestream boun- 
daries of the wake flow. To first order and for an essentially straight tra- 
jectory, p is uniform across the entire wake. To the same order, from the 
continuity equation (45), the derivative of u 2 at each freestream boundary 
is equal to the downstream derivative of Uj. Therefore, solution for 
u 2 from its momentum equation (46), as a boundary value problem is required. 

Therefore, for the cited restrictions, the developed 3DPNS differential 
equations (45)- (46) provide the required descriptions for the viscous aero- 
dynamic flows of interest. It remains to close the system with a model for 
the Reynolds stress distribution. 


Turbulence Closure and Modeling 


The operation of time-averaging has introduced the Reynolds stress tensor, 
, equation (43) into the 3DPNS equations. Using well established pro- 
cedures, e.g., reference 17, the partial differential equation system govern- 
ing the behavior of the Reynolds stress tensor can be derived. Assuming the 
cross-correlation of density fluctuations can be neglected, the exact partial 
differential equation governing the kinematic Reynolds stress tensor, -u"u* » 
is (ref. 19) 1 ^ 


L < u ?j} 


3^Vi U j } + 


— 3D i 

u j u k 3x k + 


3D j 


Vk 3^ 


„ p 

+ 2V -r r — — + — 

3X|^3Xj^ p 
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(48) 
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3x k 
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= 0 


An additional differential equation for turbulence dissipation rate, e, is 
required; under the assumption the process is isotropic, 


2v 


3x k 3x k 



(49) 
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The transport equation for dissipation function e is (ref. 18) 


L(e) = + 2v 
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3u;j 3u^ 3u£ 


3x k 8 x £ 3x £ 
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Equations (48 ) - ( 50 ) represent seven additional partial differential 
equations describing turbulence transport phenomena. This system also is 
not closed, since the third order fluctuating velocity and pressure-velocity 
correlations remain undefined. While additional differential equations could 
be established, they in turn would involve fourth order terms. Hence, a 
modeling of third order correlations is typically invoked to close the pre- 
sented system. This can be established at several levels of completeness, 
dependent upon dimensionality and geometrical complexity of the physical 
system. 

Laund er et.al., (ref. 19 ) develop a closure in terms of all components 
of -ujuj . They present results for several flow cases including isotropic 
turbulence, free shear flows, elementary duct flows and flat plate boundary 
layers. In earlier work, Hanjalic 1 and Launder (ref. 20) establish ed a 
closure applicable to thin shear flows wherein one shear component, -UjUj 
of the Reynolds stress tensor is retained, and solved in combination 
with e and the turbulence kinetic energy, k, defined as 

k = (51) 


The differential equation for k is established by contracting equation 
(48) with the Kronecker delta. For the shear dominated flows of interest, 
wherein u 1 »u 2 and u 1 »u 3 . 


L(k) 




an. 


- c 


'rz ax 


k ax, 


= o 


k 2 e" 1 


ak' 


ax, 


z 


(52) 


Equation (52) introduces the summation index convention 1 < i , j <3 and 
2 <k,&<3. The corresponding form for dissipation function, equation (40) is 
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( 53 ) 


" C e 3x i [ ke 1 U i U j 3Xj] 

9U-, _ . 

“ C ^u ^ek-^+G^k -^ 0 

In equations (52)-(53), the various constants C a are determined from approx- 
imate analyses and/or computer optimization (ref. 19)* 

The next level of simplification for turbulence closure involves specifi- 
cation of an effective kinematic turbulent diffusion coefficient, v t . From 

first principles (e.g., ref. 21), the form must be 

v t = C V Z (54) 


where C is a constant, V is a scale velocity and l is a length. For the 
turbulence kinetic energy-dissipation function, two equation closure hypo- 
thesis (herein named TKE), V is taken as the square root of turbulence kinetic 
energy, equation (51). A dissipation length scale, 2,^ , can be defined in 
terms of k and e as (e.g., ref. 21 ). 

A d = k^E' 1 (55) 

The TKE closure then yields 

v t = C v k 2 e -1 (56) 


Note that this is precisely the diffusion coefficient for turbulent kinetic 
energy, equation (52). Furthermore, upon summing the diffusion terms in 
equation (53), using equation (51) and assuming isotropy, equation (56) also 
yields the diffusion coefficient for dissipation function. 

To close the TKE system, it is required to model the correlation between 
the shear components of the Reynold's stress tensor and k, e and the mean 
velocity field u-j . For two- and three-dimensional parabolic flows, the rela- 
tion is typically of the form 


■“1“ £ s c u kV> 
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3x, 


9x i 


(57) 
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The subscript bar indicates the index is not eligible for summation. Under 
the assumption of isotropy, the elements of the correlation tensor C-j^ can 
be combined into C v , equation (56). Fjor the unidirectional flows of 
interest, partial derivatives on xi dominate. Therefore, using equations 
(56) -(57), the final form of the two-equation TKE closure is 


L(k) 


L(e) 
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(58) 


(59) 


In equations (58) -(59), for 3DPNS, the tensor indices range 1 < i , j <3 and 

2 < £ < 3. For 2DPNS , correspondingly, 1 < i,j < 2, £ = 2 only. Hence, since 
i = 1 corresponds to the direction of predominant flow, diffusion is restricted 
to the plane transverse to the Xi coordinate, as required by parabolic assump- 
tion (2). The various coefficients ( Cj<, C e , Pr G ,etc) have been determined 
from simplified analyses or computer optimization; recommended values for 
shear layer flows are given in Table 1 (e.g., ref. 19, 21 ). 

Table 1 

Coefficients in TKE Closure Model 


Variable 

Equation No. 

Coefficients 

v t 

(56) 

C v * 0.09 

k 

(58) 

Pr k = 1.0 

£ 

(59) 

Pr = 1.3, C‘ = 1.44, C* = 1.92 

£ £ £ 
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Boundary condition statements are required for the TKE equations, since 
they also represent initial -boundary value problems. In the freestream, the 
levels of k and e must vanish by assumption of a non-turbulent free stream 
flow. Since validity of the Reynolds stress hypothesis is limited to regions 
where the turbulent Reynolds number is large, it is not feasible to enforce k 
and e to vanish at a solid surface. The alternative is to express boundary 
conditions at some small distance from the aerodynamic surface. To accom- 
plish this, Prandtl mixing length theory (MLT) can be employed adjacent to 
the aerodynamic surface to close the problem, in concert with the Van Driest 
damping function distribution (ref. 17). According to equation (54), the 
effective turbulent viscosity is determined by a velocity and length scale. 

MLT expresses the correlation in terms of the predominant mean flow gradient 
and a length scale % (ref. 6 ). 


Vj. = 0) 2 & 2 


3u. 


3x 0 
L 2 J 


(60) 


where l is the mixing length 
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The Van Driest damping function is defined as 


a) = 1 - expf^A’ 1 ) 


(62) 


In equation (61), X 2 is the coordinate normal to the surface, 6 is the 
boundary layer thickness, and A and k are constants (typically 0.09 and 
0.435 respectively). In equation (62), A is a complex function of many 
factors influencing flow phenomena near the surface, including axial pressure 
gradient and normal mass flow. The form of Cebeci and Smith (ref. 17) serves 
to unify the many formulations as 
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All variables are time-averaged steady components, subscripts e and w 
refer to freestream and wall values respectively, A + is a constant (25.3), 
t w is the skin friction. Pressure gradient and mass addition effects are 
accounted for accordingly as 


v u, ) du, 
e le le 

(65) 

u 3 dx i 

U T J 1 

u" 1 V 

T W 

(66) 


where u^ is the freestream axial velocity, V w is the specified velocity 
normal to the airfoil, if present, and u T is the shear velocity 


u 


T 



(67) 


The shear stress t w is defined as 
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3D 

p w v w 3x 
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w 


( 68 ) 


It can be evaluated using a difference formulae from an empirical data-based 
equation, for example Patankar and Spalding (ref. 221, or for boundary layer- 
type flows, from the Ludwieg-Tillman formula (ref. 6~). The latter form is 


t w = |-p e u^ e [o.246(10)' 0 - 678H Re 9 ' 0 - 268 ] (69) 

where Re 0 is the Reynolds number based upon boundary layer momentum thickness 
and H is°the boundary layer shape factor (e.g., ref. 6 ). 

Equations (60)- (69) provide the formalisms necessary to determine boundary 
condition values for k and e near the aerodynamic surface. Furthermore, 
through the dual definitions of turbulent effective viscosity, equations (56) 
and (60), and since the latter involves functions only of the mean axial 
velocity component ui , which is either known or readily initialized, a means 
is established to initialize distributions of both k and e at the node 
points of a discretization. Since the developed partial differential equations 
are initial value, this information is required to start a solution. 
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The closure for turbulence is now complete at the level of sophistication 
selected for these studies. The solutions of equations (58)- (59) provide 
closure for the 3DPNS equations (41)-(42) via equation (57). Define an 
effective diffusion coefficient as 


^ E ^ + p y t (70) 

Then, using equations (70), (54) and (44), and neglecting fluid dilitation, 
the total stress tensor contribution for 3DPNS can be written as 
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Then, the 3DPNS system for isoenergetic flows equations (45)- (46) can be 
compactly written as 
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(72) 

(73) 


where 1 < i , j <3 and 2 < £ < 3. 

The viscous flowfield equation system is now complete as equations (58), 
(59), (72) and (73). This system must be solved in concert with the definition 
of the potential flow as obtained by solution of either equations (21), (25), 
(30) or (32). All boundary conditions have been identified for the repre- 
sentative solution domains. Appropriate ini tial -value character has been noted 
and means proposed to initialize required distributions in terms of readily 
available information or data. The system thus developed appears complete. 

It now remains to establish the numerical solution algorithm for these equation 
systems. 


FINITE ELEMENT SOLUTION ALGORITHM 


The desired forms of the partial differential equation systems governing 
subsonic aerodynamic flows of interest have been presented. Each member of 
the system is a special case of the general, second-order non-linear elliptic 
boundary value partial differential equation 
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Herein, q is the generalized dependent variable, the tensor indices range 
2 £ k, £ < 3 and 1 < i < 3, K is the diffusion coefficient, f^ is a function 
of its argument that specifically includes three-dimensional convection, p 
is a generalized solution parameter, and f 2 is the initial -value operator. 

The boundary condition statements for each of the dependent variables can 
be concisely expressed in the form 

£(q) = a^q + a^K^-n, + = 0 (75) 

dx i 

i.e., the normal derivative of q is constrained by q and a parameter as 
determined by specification of the a (i). An initial condition is required 

for q identified with each dependent variable as, 


q(x 1 (0),x 2 ,x 3 ) = q 0 (x 2 »x 3 ) 


(76) 


The finite element solution algorithm is based upon the assumption that 
L(q) is uniformly parabolic within a bounded open domain ft ; that is, the 
lead term in equation (74) is uniformly elliptic within its domain R, with 
closure 8R, where 


a e 



(77) 


and x 0 < X- For the 3DPNS equations, x is associated with the x^ coordi- 
nate. Equation (75) expresses functional constraints on the closure of 
ft, 3ft = 9R x [xo»x) » and the ini tial -condi tion specification, equation (76), 
lies on RU3R x Xo* 

The concept of the finite element algorithm involves the assumption that 
each dependent variable is separable in the form 


q*(x i ) = q 1 (x 1 )q 2 (x 2 ,x 3 ) 


(78) 


The functional dependence in q2^ x 2, x 3^ is represented by a polynomial 

in x^ . The expansion coefficients q^ can be most conveniently expressed 
in terms of the value of q*(x-;) at the nodes of the finite element discreti- 
zation of R. Then, equation (78) takes the form 



^l^ x £^l^ x l^ ^2 (x^, ( X x ) ^3^ x ^,^3^ x i^ 
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where the polynomials <J>(xo) are known functions of X£ and X 3 « Since they 
are known, they can be differentiated analytically, e.g., 

3q* 

= (80) 

Hence, there is no need to establish difference formulae to approximate the 
differentiated terms in L(q). 

The finite element solution algorithm is established for the equation 
system (74)- (76) using the method of weighted residuals (MWR) formulated on 
a local basis. Since equation (74) is valid throughout ft , it is valid 
within disjoint interior subdomains described by (x-j.xJeRm [xo x)» 

called finite elements, wherein UR m = R. The approximate solution for q 
within R fn x[x 0 x)> called q*(x.j , x)> 1S given in equation (79) • Therein # 
the functionals 4 >k( x &) are subsets of a function set that is complete on 
R . The expansion coefficients (x) represent the unknown x - dependent 
values of q* (x.j,x) at specific locations interior to R m and on the closure 
SRj^ , called nodes of the finite element discretization of R. 

To establish the values taken by these expansion coefficients, require 
the local error in the approximate solution to both the differential equation 
L(qfr) and the boundary condition statement &(qffi) tor 3R m n3R 4 = 0 > be 

rendered orthogonal to the space of the approximation functions. Employing 
an algebraic multiplier X , the resultant equation sets can be combined as 


{<j>(x £ )} L(q*)dr - X 
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(81) 


where S is the mapping function from the finite element subspace R m to the 

global domain R , commonly termed the assembly operator. The number of 
equations (81) prior to assembly is identical with the number of node points 
of the finite element 


Equation (81) forms the basic operation of the finite element solution 
algorithm and of the COMOC computer program to be described. The lead term 
can be rearranged, and X determined by means of a Green-Gauss theorem: 



For BRflSRm nonvanishing in equation (82), the corresponding segment of the 
closed-surface integral will cancel the boundary condition contribution 
equation (81) by identifying x a ( 2 ) with K, equation (74). The contribu- 
tions to the closed-surface integral, equation (82) whereat = 0, can 

also be made to vanish. The globally assembled finite-element solution algo- 
rithm for the representative partial differential equation system then 
becomes 
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(83) 


The rank of the global equation system (83) is identical with the total 
number of node points on RUaR for which the dependent variable requires 
solution. For f 2 non-vanishing, equation (83) is a first-order, ordinary 
differential equation. Solution of this system is obtained by COMOC using 
a predictor-corrector finite-difference numerical integration algorithm. 
For f 2 vanishing identically, equation (83) is large order algebraic, and 
the matrix structure is sparse and banded. Solution of this system is 
obtained in COMOC using a banded Stoolesky equation solver. Solution is 
also required for the continuity equation (72), which is retained for 
boundary-layer flows. Since it exists in standard form as an ordinary 
differential equation, direct numerical integration yields the required 
solution at node points of the discretization. 
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COMOC COMPUTER PROGRAM 


The finite element solution algorithm is utilized, as observed in the 
previous section, to cast the original initial -boundary value problem des- 
criptions into large-order systems of purely initial -value or algebraic 
equations. The COMOC computer program system is being developed to transmit 
the rapid theoretical progress in finite element solution methodology into a 
viable numerical solution capability. COMOC integrates or equation solves 
the discretized equivalent of the governing equation systems. Initial dis- 
tributions of all dependent variables may be appropriately specified or 
computed, and boundary constraints for each dependent variable can be spec- 
ified on arbitrarily disjoint segments of the solution domain closure. The 
solutions for each dependent variable, and all computed parameters, are es- 
tablished at node points lying on a specifiable nonregular computational 
lattice, formed by plane tri angulation of the elliptic portion of the solu- 
tion domain fi, i.e., RU9R. 

Detailed discussion on the functional design of COMOC is presented else- 
where (ref. 23). The COMOC system is built upon the macrostructure illus- 
trated in Figure 2. The main executive routine allocates core, using a 
variable dimensioning scheme, based upon the total degrees of freedom of the 
global problem statement. The size of the largest problem that can be solved 
is thus limited (only) by the available core of the computer in use. The 
precise mix between dependent variables and parameters, and fineness of the 
discretization, is user-specifiable and widely variable. The Input module 
serves its standard function for all arrays of dependent variables, para- 
meters, and geometric coordinates. The Discretization module forms the 
finite-element discretization of the elliptic solution domain and evaluates 
all required finite-element nonstandard matrices and standard-matrix multi- 
pliers. The initialization module computes the remaining initial para- 
metric data required to start the solution. The Integration module consti- 
tutes the primary execution sequence of problem solution, and utilizes a 
highly stable, predictor-corrector integration algorithm for the column 
vector of unknowns of the solution. Calls to auxiliary routines for para- 
meter evaluation (effective viscosity, Prandtl number, source terms, etc.) 
as specified functions of dependent and/or independent variables, as well as 
calls for equation solving algebraic systems, are governed by the Integration 
module. The Output module is similarly addressed from the integration sequence 
and serves its standard function via a highly automated array display algo- 
rithm. COMOC can execute distinct problems in sequence, and contains an auto- 
matic restart capability to continue solutions. 
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Figure 2. COMOC Macro-structure. 







A key feature of the developing COMOC system is automated discretization 
and subsequent refinement using the concept of a finite element "super element." 
A geometrically complex two-(and potentially three-) dimensional solution do- 
main is resolved into a number of sub-domains, the number and shape of which 
is rather arbitrary. The main purpose is to describe the geometry of the 
complete problem as the union of many elementary geometries, each of which can 
be adequately described by a few (six or eight) geometric coordinate pairs. 

This is shown in Figure 3 for a three- and four-sided super-element shape. 

These coordinate pairs are then assumed to describe quadratic surface curvature 
(hence, no inflection point allowed), the basis of which is employed to deter- 
mine a new curvilinear coordinate system wherein the super element is a regular 
straight-sided polygon. This transformed shape is then divided into equal size 
sub-domains according to an integer input by the user. The coordinate pairs 
of the sub-domain discretization are then transformed to the physical plane 
via the inverse transformation. The program automatically determines inter- 
element connectivities, and numbers the resultant node points to minimize 
matrix bandwidth. The net result is automated generation of a finite element 
computational discretization in the physical plane. The discretization can 
be made non-uniform according to a smooth progression by placement of non- 
vertex coordinate pairs of the super element description nearer the desired 
vertex node. Figure 4 illustrates the super element and resultant finite 
element discretizations for an isolated airfoil of 12 elements and 648 finite 
elements respectively. Refinement of the discretization is accomplished at 
the super element input level by merely requesting a larger number of finite 
elements to be generated in the transform plane. It must be mentioned that 
the procedure is not absolutely foolproof since the inverse coordinate trans- 
formation is not unique. In particular, placement of non-vertex super ele- 
ment coordinate pairs too near a vertex pair can induce an inverse mapping 
that is non-planar, i.e., the generated discretization overlaps upon itself. 

A computer plot of the generated discretization can be of considerable use- 
fulness in this regard. A detailed discussion on the theoretical aspects of 
this feature is included in the Appendix. 



Fig. 3. Three and Four-Sided 



Super Element Domains 
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NUMERICAL RESULTS 


The goal of the numerical evaluations is to verify finite element 
solution methodology applied to the components of a weak interaction solution 
of subsonic two-dimensional aerodynamic flow over an isolated airfoil. Since 
the finite element algorithm, equation (81), is fundamentally independent of 
the dimension of the problem space, lower dimensional results can usefully 
verify inherent capabilities and limitations. The numerical test cases, the 
results of which are presented and discussed herein, evaluate solution 
accuracy, convergence with discretization refinement, and the use of non- 
uniform finite element discretizations to enhance both. Accuracy of potential 
flow about both a NACA 4-digit and a Joukowski airfoil is documented, 
including computational application of the Kutta trailing edge condition. 
Solution of turbulent and non-equilibrium boundary layer flows, including 
adverse pressure gradient, documents computational accuracy as well as use of 
first and second order turbulence closure. Results for wake flow solutions 
document merging of dissimilar viscous and turbulent flows downstream of a 
sharp trailing edge. The final section indicates iteration among these 
various flows for a representative weak-iteration solution. 


Aerodynamic Potential Flow 


Two accuracy measures (norms) are required to be evaluated. From the 
mathematics viewpoint, for a linear, elliptic boundary value problem, such as 
equations (28) and (32), the accuracy of the numerical solution is intrinsicly 
measured in an energy norm, E( , ) 
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Since the finite element mesh is non-overlapping on R, the second form of 
equation (84) demonstrates the integrals are performed on R m and summed over 
all M. Further note that equation (84) indeed measures kinetic energy per 
unit mass, see equation (24). Denoting h as a presentative finite element 
dimension, the accuracy of the finite element solution will converge to the 
analytical solution as (cf., ref. 24). 


E(<f>-<f>*,4>-<|>*) < C^-max|<r| 2 (85) 

In equation (85), $ is the exact solution, 4>" is its extremum second 

derivative, C is a constant, and p is an integer dependent upon the 
degree of the finite element approximation function, {^(x^)} , equation 

(79). The coefficient a is a measure of discretization skewness, and can 
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be interpreted as the extremum ratio of radii of inscribed and circumscribed 
circles about each (triangular) finite element. Numerical solutions allow 
determination of the actual value of p. If close to its theoretical value, 
one can estimate with accuracy the improvement in the solution that will 
result from refinement of the discretization. 

The second and more practically important accuracy measure is pressure 
coefficient Cp , equations (33)-(36), since its distribution drives the 
weak-interaction solution. The finite element solution computes 4 >(s) directly, 
since a(3) is non-vanishing on the airfoil surface, see equation (83). The 
fundamental evaluation is then accuracy of the computed tangential derivative 
of <f> for a curvilinear coordinate s, knowing n discrete values of <j> 
at non-uniform displacements As. This is a classic problem in interpolation 
theory, and cubic splines as well as difference methods were evaluated. For 
the former, 3<f>/3s is determined at the finite element gridpoints on the 
airfoil by interpolation of f(s) assuming a spline approximation. The in- 
crements As are evaluated using Simpson's rule to integrate the differential 
equation for arc length, 
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( 86 ) 


In turn, the tangent slope required in equation ( 86 ) is obtained from a cubic 
spline interpolation of the discrete coordinate pairs defining the airfoil 
surface. 

The first requirement is to confirm order of accuracy, i.e., the exponent 
p in equation (85), and to evaluate absolute accuracy of the finite element 
prediction of f(s). The only factor driving the solution for 4 , equation 
(32), is the gradient boundary condition enforcing flow tangency, equation 
(28). This in turn is a function only of the resolution of the surface normal 
distribution, h.j * and the manner in which the finite element algorithm ex- 
tremizes solution energy E, equation (84), subject to this type constraint. 

To maintain simplicity, all computations were performed using triangular, two- 
dimensional finite elements spanned by simplex functionals. Hence, the the- 
oretical value for convergence rate p, equation (85), is two. The problem 
selected is flow over a slender parabolic-arc airfoil with sharp leading and 
trailing edges. The airfoil is simulated as a symmetric sine wave with thick- 
ness ratio ex " 1 , where X is the period and e the amplitude of the 
contour. Since the problem possesses quadrant symmetry, only a quarter- 
domain numerical solution is required. Figure 5 displays the boundary con- 
dition constraints for this domain, the basic "super-element" discretization, 
and a representative uniform finite element discretization as generated by 
COMOC. For the slender body approximation, wherein the zero normal -wash 
boundary constraint can be applied on the mean chord (which is the xi axis), 
an analytical solution exists as (eg., ref. 25 ) 
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4>(x^) = e cos 27rx 1 A*’ 1 exp(-27rx 2 X -1 ) 


(87) 


Therefore, for this simplification, both accuracy and convergence tan be 
measured exactly. 

Convergence for COMOC-generated solutions has been determined for three- 
levels of discretization refinement for 32<M<_512, where M is the number 
of triangular finite elements spanning R. Shown in Figure 6 is the plot of 
finite element solution energy E, equation (84), as a function of discreti- 
zation. Shown for comparison is a line denoting exact quadratic convergence, 
p = 2, equation (85). As observed, the convergence rate in energy closely 
approximates quadratic for the refined grid. Interestingly (cf., ref. 24), 
the convergence is from above, i.e., 


E(cf>-cf>*,<f>-<j>*) > 0 (88) 


Hence, for perturbation potential function, solution energy appears extre- 
mized by the finite element procedure according to a complementary energy 
principle. 

Solution accuracy responds favorably to the use of non-uniform discre- 
tizations. Shown also in Figure 6 is a computed solution energy for the fine 
grid case, obtained with a discretization scaled as non-uniform to enhance 
resolution in regions where gradients in <j> were largest. As observed, a 
non-uniform discretization can provide a given level of solution accuracy 
(in E) using fewer finite elements. Conversely, for a given number of ele- 
ments, greater accuracy can be acheived using non-uniform grids. This can 
be of particular economic importance for flowfield determination about actual 
airfoil contours. 

The accuracy of the algorithm for determination of Cp , was evaluated 
for a parabolic arc airfoil of chord, C = it. Shown in Table 2, as a function 

of distance from the leading edge, is a comparison between the analytical 

"" “ ' ‘ ‘ 



Two cases are shown for each discretization, corresponding to spline and 
(second-order accurate) finite difference evaluation of 3<j>/3s, respectively. 
The maximum error is less than 1 % throughout. Interestingly, the finite 
difference evaluations evidence convergence with discretization refinement 
(maximum error decreases from 0.6 to 0 . 2 %), while the spline results are 
essentially unaffected. Shown in Table 3 is a detailed comparison between 
exact and computed values of each of the three terms comprising Cp s see 
equation (36) for the coarser grid finite element solution. The major error 
is observed to be evaluation of the dominant second term. Note that the third 
term, usually discarded for linearized theory, contributes up to about 8% of 
the total Cp at mid-chord. Hence, there appears no apriori reason to drop 
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TABLE 2 


Pressure Coefficient Distribution for 
Parabolic Arc Airfoil, e/A * 0.025 


— - - - - 

Coordinate 

x/C 



Pressure Coefficient - C p 



Numerical 

, n - 13 

Numerical , 

n = 25 

Analy ti cal 

Spl ine 

F.D. 

Spl ine 

F.D. 

0.0 

0.0417 

-.0247 

-.0657 

-.0652 


-.0651 

-.0657 

0.0833 

-.1060 

-.1049 

-.1062 

-.1049 

-.1060 

0.125 

-.1449 

-.1434 

-.1453 

-.1434 

-.1449 

0.1667 

-.1818 

-.1799 

-.1823 

-.1799 

-.1818 

0.2083 

-.2159 

-.2139 

-.2167 

-.2138 

-.2160 

0.25 

-.2468 

-.2447 

-.2478 

-.2447 

-.2470 

0.2917 

-.2739 

-.2719 

-.2752 

-.2719 

-.2741 

0.3333 

-.2967 

-.2950 

-.2982 

-.2950 

-.2970 

0.375 

-.3149 

-.3135 

-.3166 

-.3135 

-.3153 

0.4167 

-.3281 

-.3272 

-.3300 

-.3272 

-.3286 

0.4583 

-.3361 

-.3356 

-.3382 

-.3357 

-.3367 

0.5 

-.3388 

-.3388 

-.3409 

-.3388 

-.3394 
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TABLE 3 


Accuracy of Individual Terms in 



The Pressure 

Coefficient 

A1 gorithm, 

Equation 

(36) 


Coordinate 

Analytical Solution 

Finite 

Element Solution 

x/C 

* 

1 

2 

3 

1 

2 

3 

0.0 

.0241 

-.0482 

-.0006 




0.0417 

.0237 

-.0874 

-.0020 

.0237 

-.0870 

-.0019 

0.0833 

.0225 

-.1245 

-.004 0 

.0225 

-.1236 

-.0039 

0.125 

.0206 

-.1590 

-.0065 

.0206 

-.1577 

-.0063 

0.1667 

.0182 

-.1906 

-.0092 

.0182 

-.1890 

-.0091 

0.2003 

.0153 

-.2189 

-.0122 

.0153 

-.2172 

-.0120 

0.25 

.0122 

-.2438 

-.0150 

.0122 

-.2421 

-.0148 

0.2917 

.0091 

-.2651 

-.0177 

.0091 

-.2634 

-.0175 

0.3333 

.0061 

-.2827 

-.0201 

.0061 

-.2812 

-.0199 

0.375 

.0036 

-.2964 

-.0220 

.0036 

-.2952 

-.0219 

0.4167 

.0017 

-.3063 

-.0235 

.0016 

-.3054 

-.0234 

0.4583 

.0004 

-.3122 

-.0244 

.0004 

-.3117 

-.0243 

0.5 

i 

.0 



-.3142 

-.0247 

.0 

-.3142 

-.0247 


* Terms 1,2,3 correspond to order of terms in equation (36) 
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the term since its evaluation is straightforward using the developed 
algorithm. 


Determination of accuracy of computed distribution of Cp was accom- 
plished by prediction of flow over a Joukowski airfoil for several thickness 
ratios and angles of attack. Initial computations utilized total potential 
function $ , commonly employed as documented in the finite element liter- 

ature. This approach was rapidly discarded, however, as numerical results 
indicated that C D solution accuracy was highly sensitive to local varia- 
tions in discretization uniformity including inclination of the element 
diagonals. This results directly from use of the vanishing gradient boundary 
condition, equation (23), over which no direct control can be, exerted. So- 
lutions using perturbation potential function have been much more successful. 
Following detailed numerical evaluations, it was determined that the "infinity" 
boundary condition, equation (27), should be placed at least one chord length 
( C) away from the airfoil, except in the leading edge region where it should 
lie 1.5c upstream. Therefore, a computational solution domain, lying 1.5G to 
2.0<C from the airfoil was selected. The airfoil contour was encircled by a 
super-element discretization employing perturbed ellipse boundaries. This 
generates a non-uniform discretization, concentrating small finite element 
domains in regions of high surface curvature and/or large field gradients, 
and provides smooth transition from fine to coarse grids. Shown in Figure 7 
is the established super element discretization, and the corresponding finite 
element discretization as generated by COMOC. The illustrated discretization 
contains 648 triangular elements. 


For pure potential flow about a symmetric airfoil at zero angle of 
attack, the solution domain represented in Figure 7 is multiply-connected. 
Numerical solution for this case verified that the finite element procedure 
will imbed an effective branch cut from the trailing edge downwind, and will 
accurately compute a planar wake trajectory. For all other cases, the Kutta 
condition must be invoked at the trailing edge and at angle of attack, a 
jump in </> occurs everywhere across the wake trajectory. Flence, the solution 
domain becomes effectively simply-connected as will always occur upon addition 
of viscous flow effects. The specification of boundary conditions reflects 
these constraints; they are also noted in Figure 7. The vanishing normal 
gradient of <f> is enforced uniformly on the "infinity" closure, while every- 
where on the airfoil surface, the normal derivative is constrained by orien- 
tation and angle of attack, see equation (28). The nodes located downstream 
of the trailing edge incl uded-angle bisector are double-numbered and a van- 
ishing normal derivative applied on both sides of the wake trajectory. The 
distributions of the (+ and -) outward pointing normals for the wake tra- 
jectory are assumed to vary quadratical ly Thus, wake curvature induces 
smooth transition from the trailing edge to angle of attack at 1.5C down- 
stream. The arbitrary constant value of $ is set (only) at the extremum 
upstream node. Therefore, circulation, lift, stagnation point and C p are 
all direct output variables and under no direct user constraint. 


Perturbation potential function distributions were computed for flow 
over a symmetric, 12% thick Joukowski airfoil at five angles of attack, 

0° < a < 8°. Angle of attack specification was accomplished through rotation 
of the freestream in the gradient boundary condition statement as discussed. 
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Fig. 7. Super Element and Finite Element Discretizations 
For Flow Over Symmetric Joukowski Airfoil 
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Computer-generated isobar distributions for the finite element solutions are 
shown in Figure 8. Note the near exact symmetry for zero angle of attack, 
even though the discretization is not symmetric, see Figure 7. This was not 
the case for use of total potential function. The solution extremum is # 
centered on the stagnation point for a = 0°and 2°. As angle of attack is 
further increased, the extremum moves into the entire lower half plane and 
becomes distributed along the trailing edge branch cut. The increase in the 
jump of <f> across the wake trajectory is also graphically evfdent as a 
increases. For all cases, note the fidelity with which the solution approx- 
imates the vanishing normal gradient around the frees tream boundary. Relative 
magnitude of perturbation velocity levels are readily discerned since equal 
increments in <f> are plotted in each case. 

The aerodynamic lift is strictly a function of circulation, r, which 
is related in an elementary manner to the jump in <f> at the trailing edge. 
Shown in Table 4 is the lift coefficient, ^computed from the finite 
element solutions of <J> 5 0°< a < 8°- Agreement with the exact 

values is good, with a nominal error of 3% essentially independent of a. 

The non-zero computed lift for zero angle of attack is one measure of the 
numerical error induced by the non -symmetric discretization. Its value is 
certainly nominal. 

The aerodynamic parameter of primary impact for the weak interaction 
solution is surface pressure distribution, from which can be determined 
pressure coefficient as well as form drag and pitching moment. Shown in 
Fig. 9 is a composite plot of C p versus a , as determined by use of the 
finite element computed distribution of tj>(s) within the C p algorithm, 
equation (36) and using the cubic spline technique. Agreement with the ana- 
lytical solution is generally good over this range of angle of attack, re- 
garding prediction of extrema, stagnation point and in the general overall 
shape of the curves (which have not been smoothed). The extremum low pressure 
on the upper surface is universally underpredicted, as one would expect, (since 
the <p solution is basically interpolator) , but the surrounding curve shape 
is good. The same holds for prediction of stagnation point. The finite 
element bound on stagnation point location, defined as the half-domain span 
surrounding the computed maximum C p , agrees well with the exact solution, 
as shown in Table 4. 

The stability of the weak-interaction solution algorithm is fundamen- 
tally a function of phenomena in the immediate vicinity of the trailing edge 
of the isolated airfoil. It is within this region that slope discontinuities 
occur in C p . Shown in Figure 10 is an expanded scale plot of computed C p 
for a = 6° including prediction on the trailing edge wake centerline. 

General agreement with data is good, although detailed differences and dis- 
crepancies are noted.. The computed values of lower and upper surface C p 

at the trailing edge are not equal. However, the curves do cross immediately 
downstream of the trailing edge. Since the $ domain is singly-connected, 
and since prediction of C p starts at the leading edge (even at angle of 
attack) as an initial-value determination, there is no mechanism present to 
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Fig. 10 Pressure Coefficient Distribution on Joukowski 
Airfoil and Trailing Edge Wake, a = 6° 
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enforce uniqueness at the trailing edge. Shown in Table 5 is the computed 
difference in C p at the trailing edge node location as a function of a . 


Table 5 


Pressure Coefficient Difference at Trailing Edge 
Of Joukowski Airfoil 


Angle of 
Attack a 

Pressure Coefficient 

, c p 

Upper 

Lower 

Difference 

0 

0.1166 

0.1516 

-0.035 

2 

0.1094 

0.1587 

-0.049 

4 

0.08915 

0.1664 

-0.077 

6 

0.1039 

0.1700 

-0.066 

8 

0.0876 

0.1680 

-0.080 


Note that, even though the finite element solution for $ was symmetric 
to = 0.001 for zero angle of attack, that accumulated error in determina- 

of C produced a trailing edge difference of 0.035. This could be attribu- 

r 

table to the difference in diagonal orientation within the discretization, 
see Figure 7, since "averaging" solutions obtained with the given discretiza- 
tion, and one with the opposite diagonal orientation, would obviously eliminate 
the Cp difference. Alternatively, it could be accumulated round-oVf error 

in the C p spline algorithm, since upper and lower surface normal vectors are 

of opposite sign and no advantage has been taken of the intrinsic problem 
symmetry. Finally, it could be associated with the Joukowski cusped trailing 
edge, which is only approximately modeled for any discretization. From Table 
5, the computed trailing edge difference appears essentially insensitive to 
angle of attack. 

The remaining difference may be attributed to the influence of the 
assumed wake trajectory curvature on the solution. An airfoil itself provides 
ample evidence that flow curvature can support a pressure gradient. The 
exact solution, obtained using conformal mapping does not 'tDossess" a wake, 
and therefore is insensitive to the manner in which the velocity vector at 
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the trailing edge is rotated to frees tream orientation. However, in the 
present boundary value solution, the wake trajectory is as much an impervious 
surface as is the airfoil itself. Since the distribution of the normal to 
the wake is assumed quadratic, the wake trajectory itself is cubic. The 
constraints on the cubic curve are known slope at each end (Kutta condition 
and freestream), known trailing edge intercept, and the assumed quadratic 
variation of ^ • Several other trailing edge orientations were evaluated 

for the 6° angle of attack case, including a straight wake trajectory at the 
Kutta angle and at the freestream angle, and a linear average of these. In 
all cases, the pressure difference computed at the trailing edge exceeded the 
present case by a factor of up to four. Therefore, the quadratic normal 
distribution appears to best simulate the experience with exact solutions. 

Considerable additional numerical experimentation could certainly be 
conducted. The developed spline C p algorithm may not be optimal, especially 
since it did not clearly demonstrate convergence under grid refinement for 
the parabolic arc airfoil test case. Refined grid potential flow solutions 
should be generated to allow measurement of convergence in the energy norm. 

This development study could be readily conducted using the COMOC computer 
program. 


One additional point warrants comment. The design trend is to pro- 
gressively thicker airfoils. As the thickness increases, the placing of flow 
tangency boundary conditions on the actual surface (instead of the mean chord) 
becomes additionally important. The finite element concept retains all the 
required flexibility. For example, shown in Figure 11 are finite element 
discretizations for solution of 24% thick Joukowski airfoil, accounting for 
and neglecting the influence of thickness. Figure 12 compares the computed 
Cp distributions, as obtained by applying identical boundary condition 
constraints, equation (28), on the surface and on the mean chord. The 
differences are most pronounced as expected in the region of maximum thickness. 
In the same manner as produced the "flat plate" Joukowski airfoil approxima- 
tion, the discretizer in COMOC can readily establish an arbitrary non- 
symmetric airfoil shape including camber. 
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Boundary Layer Flow 

A numerical test program was employed to ascertain the accuracy and 
convergence character of finite element prediction of turbulent two-dimen- 
sional boundary layer flows using the 2DBL option in the COMOC computer 
program. Emphasis was placed on assessing the influence on accuracy and 
computer cost of finite element discretization and refinement, turbulence 
closure modeling, and initialization procedures for required data profiles. 
The comparison bases were computed mean flow longitudinal velocity (5'i) 
profiles, the familiar boundary layer integral parameters of displacement 
thickness (6*) and momentum thickness (0), and the point norms of shape 
factor (H), boundary layer thickness (<$), and skin friction (Cf). For 
reference, these definitions are 


6 * = 


J 

6 


[1 - Ve 1 ]^ 



6 

H e 6*e-i 

6 = {x 2 |Vx 2 , u^u” 1 = 0.991} 



(87) 

( 88 ) 

(89) 

(90) 

(91) 


where x 2 is the coordinate normal to the aerodynamic surface, u e is the local 
freestream velocity, and t w is the wall shear stress. Each of the defined 
parameters is a function of Xj, the longitudinal (streamwise) coordinate. 
Furthermore, for three-dimensional parabolic or boundary layer flows, each 
also depends upon the lateral coordinate x 3 . The Ludweig-Ti liman formula 
(ref. 6 ) for is valid for laminar, transitional and turbulent flows, and 
is determined from 0 and H as 


C 


f 


0.246(10)~°' 678H Re '°' 268 

0 


(92) 


In equation ('92 ), Re^ is the momentum thickness Reynolds number 


Re e 



( 93 ) 
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Four test cases were selected to facilitate the study. The turbulent 
flow test cases were taken from data sets published in the Stanford conference 
proceedings (ref. 26 ). The Wieghardt flat plate flow (IDENT 1400) corresponds 
to equilibrium turbulent boundary layer development at zero pressure gradient. 
Constant freestream velocity (U ro = u e ) was 33 m/s, the flow was tripped at 
the blunt leading edge by a trip wire, wind tunnel turbulence level was 0.25%, 
and the average unit Reynolds number was 2.19xl0 6 per meter. The Bradshaw 
relaxing flow (IDENT 2400) corresponds to evolution of a non-equilibrium 
boundary layer, induced by abrupt removal of a moderate adverse pressure 
gradient from an initially equilibrium flow. Nominal freestream velocity 
(Uoo) was 33.5 m/s, wind tunnel turbulence was less than 0.1%, and the reference 
unit Reynolds number was 2.38xl0 7 m . The Newman adverse pressure gradient 
flow (IDENT 3500) corresponds to turbulent boundary layer flow approaching 
separation on the upper surface of an airfoil at angle of attack. Nominal 
reference velocity (Uoo) was 36 m/s, wind tunnel turbulence level was approx- 
imately 1%, and the reference unit Reynolds number was 2.5xl0 6 m" 1 - Finally, 
the zero pressure gradient, laminar incompressible (Blasius) similarity 
solution (ref. 6 ) was employed to assess absolute accuracy and validate 
numerical convergence for the finite element discretizations determined econom- 
ically accurate for turbulent flow prediction. 


The two-dimensional, boundary layer equation system solved by the 2DBL 
option in C0M0C is established in a transformed coordinate system. Coordi- 
nate stretching in the x 2 direction can efficiently compensate for bounday 
layer thickness growth. Referring to Figure 13 , the selected transformation 
is 

Xp - f i (x., ) 

n = (f 2 Cx 1 ) - f 1 (x 1 ))/f (94) 


where n is the normalized coordinate lying between (at least) piecewise con- 
tinuous surfaces, fi(xi) and f 2 ( x i) > f x , that bound the elliptic solution 
domain R, see equation (77), and f is a normalizing factor. For the presented 
results, fj corresponds to the plate surface and is a constant, while f 2 
is appropriately specified for each case. (Note that for the Blasius case, 

f 2 (x ! ) ~ C\/x| would exactly scale the finite element-spanned domain to that 
of the similarity solution. However in C0M0C, f 2 is retained in complete 
generality as piecewise linear interpolation between input data points, or 
as an input functional.) 

The boundary layer equations are a sub-set of the parabol ic^Navier- 
Stokes system , equations (72)-(73). For two-dimensional flows, U 2 is de- 
termined from equation (72) and Ua is solved from equation (73) for i = 1, 

1 < j < 2 and l = 2. Equations (58) -(59) are solved for the TKE closure 
modeling (if used) over the same index range. In expanded form in the 
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(£,n) coordinate system, and applying the incompressible boundary layer order 
of magnitude simplification, the equations solved are 


L(u 2 ) 

Map 


= [35 ■ ^ h 2 + nh 3 ) 8n] a i 
= u l[w " ^ h 2 + nh 3^ u l 


au 9 

+ -W-° 


+ u 


3iL 

2 IrT 


(95) 

(96) 


_3_ e. 

3q 1 3q 


1 d P( 
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L(k) = 


p 

D i[w' (h 2 + nh 3 ) lr] l< + a 

(v e h 


3k 
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_3_ 
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1 3k 
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e K] 2 
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L(e) = u- 


M " (h 2 + 1 h 3hsr 


8n 


. ~ 3e 
e u 2 3n 


(98) 
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3n 


v e h 


1 3e 


Pr e 3q 



3u n 
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1 

. 39, 


+ C 2 e 2 k” 1 = 0 
e 


The functions h^ , 1 ^ i ^3, are the metric coefficients which result from 
differentiation in the curvilinear coordinate system defined in equation (94) 
Members of the set are defined as f 

m 

2 


f(f, - V 1 


h. h {f'h 


1"1 

h l (f 2 ' 


> 


(99) 


The superscript prime denotes ordinary differentiation with respect to 
The effective kinematic diffusion coefficient is (see equation 70). 


v e = Re” + v 


(100) 
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where 




3iL 


v t =< 

(i 

1 

1 

3x 2 

(MLT) 


C v k 2 e _1 

(TKE) 


(101) 


and v is the kinematic molecular viscosity. For the MLT closure, Z is the 
mixing length and a) is the Van Driest damping coefficient (ref. 17), see 
equations (61) -(62,). 


As the key feature of turbulent flow computations, it is mandatory to 
employ a highly non-uniform (finite element) discretization of R to obtain 
satisfactory computational efficiency in concert with acceptable solution 
accuracy. Following extensive numerical tests, it was determined that 
solution speed and accuracy were enhanced in concert by using a finite element 
discretization that increased in sequential span according to a geometric 
progression. For two-dimensional parabolic flows, the finite element domain 
Rro is a line segment of measure (length) L m . If n m= i represents the ex- 
tremum nodal coordinate of i.e., the mill finite element domain spans 


L m E ^m+l - 


( 102 ) 


then a geometric progression of n (hence span of Rj^) is obtained by the 
recursion relation 


n m+l = n l + s .I p ° l<m<M (103) 

j 

In equation (103), ru is the coordinate of the first node of the discre- 
tization (typically zero, i.e., x 2 = fi(xi) , equation (94)) p is the pro- 
gression ratio of the finite element span, s is a scale factor that allows 
imbedding a given number of elements within fz - fi, equation (94) , and 
M is the total number of finite elements spanning R, equation (77). 

Numerical experimentation was conducted for 

10-3 < s/6 < HT 1 
1 < p < 1.5 

where 6 is the boundary layer thickness. Note that p = 1 corresponds to a 
uniform discretization. Shown in Figure 14 are graphs of discretizations 
as established by equation (103 ) for several s/6 and p. Curves A, C and D 
illustrate various uniform discretizations of R; the header shows the corres- 
ponding number of finite elements M spanning 6 and 36, whereat appropriate 
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freestream boundary conditions for each computed variable were applied. Curve 
B illustrates a modestly non-uniform grid, suitable for laminar flow predic- 
tions, obtained for p - 1.11. Compared to curve A, resolution at the wall 
is doubled, while only about two-thirds the number of elements is required to 
span the solution domain, R = 36. Curve E depicts a typical finite element 
discretization determined by numerical experiment to yield good solution 
accuracy for turbulent flow predictions while requiring only minimal computer 
execution time. The finite element at the wall spans 6 x 10 " 3 , yielding 
excellent resolution of near wall damping, yet only about 28 finite elements 
are needed to span R = 36. Progression ratios p > 1,25 generally yielded a 
too sparse discretization in the freestream. Values of p appreciably less 
than 1.25 can place double the number of elements within 36. Barring the 
numerical round-off error associated with finite arithmetic, finer discreti- 
zation always improves accuracy for a consistent algorithm. However, for the 
turbulent flows tested, the improvement in accuracy was typically inadequate 
to compensate for the doubling of computer CPU time. Hence, the following 
results were obtained using the discretization shown as curve E, Figure 14 

The results for the Wieghardt data case (IDENT 1400)_are presented in 
Figures 15-16. As shown in Figure 15 the computed Ui velocity profiles 
are in excellent agreement with data at each of the first eight experimental 
stations. Closure for turbulence was accomplished using MLT, with the para- 
meters k and X equated to their standard values of 0.435 and 0.09 respectively. 
For boundary conditions, both ui and U 2 vanished at the plate surface, while 
3Ui/9x 2 vanished at the freestream boundary, x 2 = 36. Shown in Figure 16 
are comparisons between data and the finite element predictions for the 
boundary layer parameters. Agreement is excellent throughout, except for skin 
friction C^ , where measureable disagreement is noted over 0 < 0.25m. 

However, within this range, the experimental data is extremely sparse; note 
the lowest symbol for ^ on each data profile in Figure 15, corresponding 
in each case to the first measured value. Since Cf is basically the wall 
shear stress t w , its evaluation is highly sensitive to the employed inter- 
polation formulae (see comments, ref. 26, Vol . II, p. 98). 

The computed solution was initialized in a manner identical to experiment 
wherein a turbulence-free uniform flow was forced to impinge upon the plate 
leading edge. For the numerical solution, the initial value of 1 ^ at each 
node of the finite element discretization was thus uniformly Uoo, ' except for 
the node at the plate surface where Uj = 0 throughout. Hence, the maximum 

attainable velocity gradient occurred at solution initiation; nevertheless, 
start-up by C0M0C was smooth and devoid of oscillation. The initial values 
for U 2 were uniformly zero since aui/ 3 xi = 0 upstream of the plate. This 
corresponds to what is termed a "slug-start," which is the method employed to 
initiate solutions in the absence of any preferable alternative. Since the 
turbulent Reynolds number of the Wieghardt case is initially small, transition 
from laminiar to turbulent flow occurs over a finite span of x 1 , up to and 
including perhaps the first data station. The results of computational ex- 
periment with the switch from laminar to turbulent flow is shown in Figure 17 
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The symbols correspond to Cf as established from Wieghardt data using both 
an interpolation procedure (ref. 26) and the kudweig-Ti liman formula, equation 
(A. 6). The curves were determined from the ffj- profiles computed by COMOC 
using equation (A. 6). In each case the flow was assumed laminar to the 
selected transition point after which computation of the MLT eddy viscosity 
was initiated and continued as the flow proceeded downstream. No transition 
model was employed to alter the intermediate profiles. The family of computed 
results are bracketed by the data, and by x x = 0.6m the various methods are 
in essential agreement. The velocity profiles shown in Figure 15 were taken 
from the run which predicted uj at xi= 0.087 m in agreement with the first 
experimental value ( g .77). 


The second turbulent flow test case is considerably more demanding* 
since nonequilibrium phenomena are involved in the relaxation process. Tur- 
bulence generation and dissipation processes are not in balance at a point in 
the flow; hence, one expects MLT to be of marginal validity while the TKE 
closure should be more exact. Nevertheless, MLT with standard coefficient 
values proved remarkably capable as illustrated by the detailed Uj velocity 
comparisons in Figure 18 These results were generated using a fixed coor- 
dinate system since boundary layer growth was rather modest for this case. 
Agreement with data is good, although diffusional processes within the flat 
mid-range appear too high as evidenced by the computed results uniformly 
exceeding the data. These differences diminish further downstream, but 
there is a corresponding trend for the computed results to under-predict the 
first knee in the curve. Shown in Figure 19 is comparison between prediction 
and data on the boundary layer parameter bases. As suggested in reference 26, 
match with the data has been enforced at the second station. Agreement in 
0 , H and 6 is generally excellent.' The correct trends and local extrema in 
6* are predicted; however, the overall level of the curve is high. In concert 
with the underprediction of Cf> these curves are in agreement with the computed 
low profiles near the wall. Note that the increasing wall shear in a 
uniform pressure field is a manifestation of the non-equilibrium processes. 
Turbulence levels for this case are high as well. For example, see equation 
(100 )> the computed extremum of ve/v for the Wieghardt case was 150, while 
for Bradshaw it was 900. However, for both cases, the employed non-uniform 
discretization provided adequate resolution of the wall region while using 
only 29 elements to define the range of extremum turbulent viscosities 
encountered. 


The Bradshaw case was employed to assess correct operation of the TKE 
closure model for turbulence in COMOC in concert with MLT to predict near 
wall effects. Several tests were run to evaluate an optimum y+ to apply 
fixed boundary conditions for solution of k and e, equations (97 )-( 98 ), 

where y+ is the familiar Reynolds number based upon u , equation (67), 
i . e . T 
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Initial profiles for k and -e , to start the solutions, were determined from 
the Uj profiles using the dual, definition of * equation (101), in concert 

with the definition of a dissipation length scale, 

l d = k^e -1 (105) 

and the knowledge that SL^ can be interpolated quadra tically between known 

wall and freestream values (0 and 0.09, respectively), and that 3£ C |/3X2 at 
the wall is linear and known (e.g., see ref. 21). The range evaluated was 
7 < y + < 50 » and. best results for this case occurred for 30 < y + < 40* 

Shown in Figure 20 are comparisons between computational results for Uj 
and data; agreement is good and comparable to the MLT run. Figure 21 pre- 
sents comparisons for the boundary layer parameters. Agreement of computed 
< 5 * and 0 with data is noticeably poorer for Xi > 71, but otherwise in essen- 
tial agreement with the trends noted for the MLT runs. 

The evaluation of the Newman adverse pressure gradient case (IDENT 3500) 
indicates limitations of the boundary layer and mixing length concepts. This 
flow corresponds to upper surface boundary layer development on an airfoil 
at angle of attack sufficient to induce separation before the trailing edge. 
Using the standard values for the constants in MLT, k = 0.435 and X ~ 0.09, 
and the previously verified non-uniform discretization, numerical evaluation 
indicated that the computed eddy viscosity was excessively large, since the 
velocity profile shape near the wall was considerably less retarded than the 
experimental data. Retardation of the boundary layer at freestream was 
correctly computed. Since the adverse pressure gradient acts as a sink for 
momentum across the entire boundary layer, i.e., it is a uniform solution- 
insensitive non-homogeniety in the differential equation, the only means to 
under-decelerate the wall profile is for the high freestream momentum to 
cascade into the low momentum region. Since u 2 «u 1 , the sole mechanism 
producing this phenomena is an excessively large diffusion coefficient, y e . 
Since some latitude exists for X , a value near its lower bound (0.07) was 
evaluated. The results of this prediction are shown in Figure 22, which 
compare computed u 1 velocity profiles to data over a one meter span. Agree- 
ment is excellent up to the last two profiles; over diffusion progressively 
erodes this agreement thereafter and the separation imminent at the last data 
station is not captured. This results directly from excessive momentum 
diffusion to the wall. For example, shown for comparison as the dashed curve 
is one velocity profile computed for X = 0.09. For this case, the extremum 
■y e /y was 820; for X = 0.07 at the same station, maximum y e /y was 650. 
However, this level increased to 850 by the last data station, resulting in 
the indicated poor agreement. The comparison of computed integral parameter 
distributions to data is shown in Figure 23. Agreement is excellent up to 
the imminent separation. Thereafter, only 6 agrees well with data, indica- 
tive of correct simulation of the pressure gradient effects using the vanish- 
ing gradient freestream U! boundary condition. This type of boundary layer 


63 



Longitudinal Velocity 


BRADSHAW 



rig. 20. Longitudinal 


FLOW (IDENT 2400 


53.0 59.1 65.1 71.1 83.0 

53. 59. 65. 71. 83. 



Fransverse Coordinate - x 2 /<$ r 


locity Profiles, Bradshaw Relaxing Flow? TKE 




Longitudinal Velocity 










situation is always encountered near the trailing edge of practical airfoil 
flows. Careful attention to detail is therefore required, including perhaps 
use of a higher-order closure model for turbulence like TKE. 

Absolute accuracy and convergence character of the finite element algo- 
rithm at this point can only be determined by computation of laminar flows. 

The Blasius solution (ref. 6 ') provides the required exact solution to 
equations (95)-(96) for v equal to a constant. Use of a nominal 10 element 
uniform discretization, curve A of Fig. 14, produces a maximum error in ui 
of about 1.5%. The non-uniform discretization of the same number of elements, 
curve B, reduces this error by a factor of about four, as does a doubling of 
the uniform discretization to 20 elements. The latter results confirm that 
the linear functional finite element algorithm is of second-order accuracy. 

It is also noted that use of the variable stretched coordinate option, 
equation (94), degrades absolute accuracy for laminar flows by up to a factor 
of two. The culprit appears to be the continuity equation solution, which 
as a direct consequence of the transformation contains a term that is linear 
in n . Hence, it grows without bound as equation (95) is solved. Use of 
the Blasius solution for u 2 , in concert with integration of equation (96) 
for u, is determined to remove the additional error. This phenomena appears 
of negligible influence for computation of turbulent boundary layers for the 
cited test cases, provided sufficient resolution of the region near free- 
stream is retained. Recall that large finite elements are therein employed; 
use of coordinate stretching can induce a too coarse discretization near 
freestream for computations of sufficient length. 


Turbulent Wake Flow 

Prediction of the merging of dissimilar turbulent boundary layer profiles, 
downstream of the terminus of the airfoil, requires solution of the 2DPNS 
equation system closed with the TKE model. The primary requirement in this 
solution is to accurately predict the rapid reduction in boundary layer dis- 
placement thickness, defined now as the sum of the upper and lower surface 
contributions, as the viscous flows merge. In an interaction solution these 
data will reflect back into the potential flow solution as alteration of the 
effective inviscid shape of the airfoil and its wake. The equilibration process 
will occur very rapidly within about 0.1.C downstream and will locally generate 
large levels of turbulence in the initial momentum— defect region. Within this 
span, the initially zero u x at the 'trailing edge will rapidly accelerate to 
50 - 60% of U ro . As the flow proceeds further downstream, the momentum defect 
within the wake will continue to slowly diminish. 

From the aerodynamics standpoint, the output of primary value for the wake 
flow is prediction of the drag due to friction. An analytical exact expression 
for drag can be derived from the integral momentum equation (cf., ref. 6). For 
an incompressible flow, the form is 
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where 0 and H are parameters defined in equations (88)- (89). In the 
absence of a wall, this equation is formally integrable by parts, yielding 
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where 0 and H may be evaluated at any x 1 coordinate. The trailing edge 
has been chosen, typically, since most interaction solutions have available 
only the boundary layer parameters ud to this location. Evaluation of 
equation (107) requires G^e = u ie (H). In the Squire-Young approximation 
(cf., ref. 27), a linear relation in the logarithm is chosen. The resultant 
drag coefficient is 
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where 0 and H are evaluated at the trailing edge. The Squire-Young 
approximation is at best only marginally valid between the trailing edge and 
downstream infinity. However, its use may be of greater value if the transi- 
tion from the trailing edge to some nearby downstream station could be more 
accurately predicted, i.e., numerically. Then equation (108) could be employed 
to complete the computation of Cq for the remainder of the wake. 


The basic requirement, therefore, is to use the 2DPNS solution to 
compute the near wake. The Joukowski airfoil test case was employed since 
experimental data base exists (ref. 28). The measured u 1 boundary layer 
profiles were interpolated to provide an initial condition for the 2DPNS 
equations. Starting at x x /C = 0.98, the 2DBL_option in C0M0C was run using 
MLT, then the TKE closure to equilibrate the u 1 and computed u ? profiles 
with the imposed C p . At Xj/C = 1.00, a restart was used to switch to the 
2DPNS option. Thi s^sol ution was then marched downstream to x^C = 1.25. 

Shown in Figure 24 is the comparison between Ui distributions and u 2 
distributions computed by the finite element 2DPNS. Agreement is excellent 
to x /C = 1.10; thereafter, the upper profile is progressi vely underdiffused. 
Note that the wake is modestly curved, as defined by the trajectory of Uj 
minimum. While the velocity differences are only a few percent, the impact 
on integral parameters is pronounced as shown in Fig. 25 for the upper 
profiles at Xj/C > l.l. The lower profile integral parameter solution is 
in good agreement with data as are both shape factor contours and the 
acceleration of the Q A minimum. Shown in Fig. 26 are computed wake dis- 
tributions of turbulence kinetic energy and dissipation. Assuming the 
validity of equation (47), the TKE profiles correspond to pressure variation 
across the boundary wake. The initial levels are modest with a double peak. 
They coalesce into a significant extremum, immediately downstream of the 
trailing edge, with higher levels maintained thereafter, coincident with a 
gradual broadening of the base. These results verify the fundamental 
capability of the finite element solution of the 2DPNS equations to accurately 
predict the merging of dissimilar turbulent boundary layer profiles within 
the wake downstream of a sharp trailing edge. The results are applicable to 
both establishment of a weak-interaction solution algorithm as well as pre- 
diction of viscous drag. 69 
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An Interaction Solution 


The finite element solution of the fundamental building blocks of a 
weak-interaction algorithm for subsonic aerodynamic flows has been documented. 
An evaluation of these capabilities for a complete flow solution has been 
completed. In this case, the airfoil contour was a modified NACA 0015, for 
which experimental data are available (ref. 29). Shown in Fig. 27 is the 


COMOC-computed C. 
tional results ar! 


distribution for this airfoil at a = 8°. The computa- 
in good agreement with the exact solution, including the 


trailing edge, where AC p - 0.01 between lower and upper surface. In dis- 
tinction to the JoukowskY, the NACA 0015 has a finite thickness trailing 
edge. Use of the original finite element solution-domain diagonal orien- 
tation, Fig. 7, for this case is thus observed to produce excellent results. 


C0M0C contains a subroutine to integrate the integral momentum form 
of the boundary layer equations. This solution is quite inexpensive and can 
yield adequately accurate determination of e and s* up to the immediate 
vicinity of the trailing edge, i.e., x 1 /C < 0.9. The upper and lower surface 
Cp distributions computed by C0M0C for the basic airfoil were employed to 

compute the initial estimate of < 5 * distributions. The trailing edge values 
Were then faired into a representative freestream decay, and the potential 
flow determination repeated for this viscous-augmented shape. This Cp 

distribution is also shown in Fig. 27, denoted "augmented '(0) ," ie., the 
zeroth iteration. Note the pronounced blunt base introduced at the airfoil 
trailing edge by the viscous effects, with an attendant sharp slope disconti- 
nuity in proceeding into the wake. 

Also shown in Fig. 27 is the 6* distribution computed from the initial 
(inviscid) computed C p and that of the first iteration. The results at the 

trailing edge differ remarkably, with the upper surface peak in 6* moving 
upstream from the airfoil terminus. The iterative cycle between viscous- 
augmented potential flow and integral boundary layer and wake was continued 
until convergence was achieved, defined as AC p $ 1% of the maximum C p . 

The final results for C p and 6* computed consistent are also shown in Fig. 

27. The latter curve is in essential agreement with experimental data (ref. 
29). 


This converged C p distribution is then employed to drive a 2DBL and 
2DPNS solution from upstream of the trailing edge on into the turbulent wake. 
Computations are continued to x^/C - 1.15, whereafter <S* and 0 can be 
assumed to vary only slowly to the end of the computational domain. The 
computed modifications to S* within the interval 0.9 < x^/C <1.15 are 
then employed to modify the effective inviscid airfoil shape. The finite 
element determination of potential flow may then be repeated if sufficient 
difference is computed for <$* from the more exact solutions. The 2DPNS 
wake flow solution is also used to evaluate drag, equation (107). 
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CONCLUDING REMARKS 


The results of the reported comprehensive numerical experimentation 
document the basic capabilities of finite element solution of the differential 
equation systems constituting two-dimensional weak-interaction theory in 
aerodynamics. Factors affecting solution accuracy for potential flow deter- 
mination about an isolated elementary airfoil have been evaluated. An algo- 
rithm for accurate computation of pressure coefficient is established. A 
highly non-uniform discretization has been determined to allow economical ly- 
accurate prediction of two-dimensional turbulent boundary layer flows with 
pressure gradient. Comparisons have been made with experimental data using 
both MLT and TKE closure for turbulence. The prediction of turbulent merging 
of dissimilar turbulent boundary layer profiles has been accomplished in the 
wake downstream of a sharp airfoil trailing edge. The ability to combine 
these analyses into a weak interaction solution for the complete flow about 
an isolated elementary airfoil has been documented. 

These results confirm for the first time the fundamental capability of 
finite element solution methodology applied to the complete aerodynamics 
problem. There is ample room for improvement as well as extention to more 
complex airfoil systems and higher dimensional space. A systematic numerical 
experimental project should be persued to evaluate under close control the 
factors affecting accuracy of computed potential function and derived pressure 
coefficient. A higher-order, curved-sided finite element might prove eco- 
nomically useful to more accurately model airfoil surface and wake curvature. 
The accuracy of the turbulent boundary layer and merged wake flow solutions 
could be enhanced through turbulence model refinement as well as use of higher 
order accurate finite element functionals and initial -value integration pro- 
cedures. Each of these areas could and should be addressed to render tangible 
the potential usefulness of the finite element solution in subsonic aerody- 
namic system design. 
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APPENDIX 


Automatic Discretization Refinement 

Combined viscous-inviscid flowfield iterative solutions impose*, the 
requirement for automatic redistribution of the potential flowfield grid 
points to accomodate the inviscid boundary shift at each solution iteration. 
The finite element size along the boundary must remain relatively constant, 
however, to maintain solution accuracy, hence, the entire flowfield grid 
must be shifted to conform to each new boundary shape. For the infinite 
variety of airfoil, boundary layer and viscous wake shapes encountered in 
aerodynamics, a method is required which is unrestricted in its ability to 
model general geometric shapes. 

The COMOC discretizer has this capability. By subdividing the flow 
domain into a few shapes having quadratic variation along the boundaries, 
virtually any closed domain can be accurately modeled and discretized to 
conform to solution requirements. The method is graphically presented in 
Figure A.l for a single element two-dimensional airfoil in an infinite 
domain. In Figure A. la the flowfield has been divided into twelve sub- 
domains consisting of eight quadrilateral and four triangular shapes. The 
gridpoints represent the data required to specify the sub-domain boundaries. 
Note that all boundaries are of degree 2 or less. The side gridpoints serve 
to describe the curved surface shapes and their location relative to adjacent 
verticies determines the spacial distribution of the generated grid. Figure 
A. lb illustrates a 648 element grid generated from the subdomain description 
in Figure A. la. Note that the refined grid around the airfoil surface and 
toward the leading and trailing edges is due to the off center placement of 
the subdomain side gridpoints. Variations on refinement type and number of 
elements is readily accomplished from a single set of subdomain data by 
simply moving side gridpoints and respecifying the number of elements for 
each subdomain. 

The grid refinement algorithm proceeds by looping over the subdomain and 
generating gridpoint and finite element correction data. For each of the 
quadrilateral and triangular subdomain shapes, a coordinate transformation- 
exists which maps a local natural coordinate system onto the physical plane. 

The general form of the required transformation for all space is: 

x n -= Qj{Xj } j = 1, n (A.l) 

where n is the total number of specified grid points on each subdomain boundary 
and summing over repeated indices is observed. The scalar x-j consists of the 
gridpoint physical values of the coordinates at the i th generated node. The 
quantity Qj is called the shape function and contains functions of e and 
n which satisfy the transformation illustrated in Figure A. 2 for e and n 
of the {xj}. Equation (A.l) is valid over one-, two-, and three-dimensional 
cartesion space and shape factors may be derived for a variety of geometric 
shapes and polynomial degrees. Since the derivation procedures are well 
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documented (ref. 8), the two-dimensional quadrilateral having a biqua- 

dratic (8 coef.) interpolation function will be focused upon with the under- 
standing that other shape functions may be directly substituted. The quadri- 
lateral shape function therefore is: 

'k{l - e)(l - n)(-e - n - 1) ' T 
h(l + e)(l - e)(e - n - 1) 
k{l + e)(l + q)(e + q - 1) 

q. = J %(.l - e) (1 + q)(-e + n - 1) l (A. 2) 

3 %(l - e *)(l - n) 

h(l + e)(l - n 2 ) 
h(l - e 2 ) (1 + n) 
h(l - e) (1 - n 2 ) 



A. 2a Physical Plane A. 2b Transformed Plane 


Fig. A. 2 Mapping of a Quadrilateral Onto a Natural Coordinate Plane (e q) 
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where the e, n axis bisects the quadrilateral and the equations are ordered 
according to the gridpoint numbering noted in Figure A. 2. Substitution of 
equation A. 2 into A.l for a specific set of {xj} and evaluating the equation 
over the limits of e and n (-1 to 1) yields a biquadratic approximation 
to x over the subdomain. Accuracy of the values of x-j are dependent upon 
the ability of the shape functions to approximate the physical geometry. A 
curvature which is exactly biquadratic, for instance, will be interpolated 
exactly using equation A. 2. Note in Figure A. 2b that the side located grid 
points lie at exactly midside, by definition. The side nodes in Figure A. 2a, 
however, need not be at exactly mid-side since it is not required in the 
definition of {Q}. Movement of the specification of the side nodes in the 
physical plane, therefore, allows for a smoothly varying distribution of 
generated data over the subdomain as illustrated in Figure A.l. The coor- 
dinate transformation is of the serendipity family and can yield non-unique 
distributions if excessive skewing of the mid-side node is attempted; there- 
fore, care must be taken when applying the method as in Figure A.l, where 
regions of highly refined grid are necessary along the airfoil surface. 

Each of the subdomains of which the entire solution domain is composed 
may be treated independently if the data generated along subdomain boundaries 
is identical for each subdomain. This requires that inter-subdomain boun- 
daries share the same gridpoints and that (e,n) be independent along the sub- 
domain interfaces. To illustrate that this is indeed the case, consider the 
two coincident subdomains of Fig. A. 3 

Along the common boundary (3,4,7)* and (4,1,8) 2 


n x = c = l.o 

= c = -1.0 

and evaluating equation (2) with these constraints yields: 

Q * L%(c + e 2 ), h{-e + e 2 ), (1 - e 2 ) (A. 3) 

Q = lh (-n + n 2 )> %(n + n 2 ), (l - n 2 ) (A. 4) 

And since -n 2 along the boundary, equation (A. 3) becomes 

Q = L%(-n + n 2 )» h{r) + o 2 ), (1 - n 2 ) (A. 5) 
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Fig. A. 3 Two adjacent subdomains 


which is identical to equation (A. 4). Other boundary combinations and the tri- 
angular shaped subdomain can be shown to yield similar results. Three- 
dimensional subdomains have surface interfaces and similar equations can be 
derived as above where the equivalent of equation (A. 5) contains functional 
expressions for a planar surface. This demonstrated generality, therefore, 
completely eliminates any requirements for specifying the subdomains in a 
particular order and numbering may begin at any vertex. 

It is impossible to encompass the entire spectrum of geometric shape and 
refinement requirements within a few subdomain types. The completeness of 
a data refinement scheme, therefore, requires that addition to its geometric 
flexibility be open ended. Utilizing the subdomain method, special subdomain 
functions may be developed and inserted directly into the system. An example 
of this is a telescoping quadrilateral having quadratic interpolation. The 
purpose of this special function is to provide a means for increasing or de- 
creasing the number of generated triangular elements when progressing from 
one end of a subdomain to another. Fig. A. 4 illustrates the effect for a 
quadrilateral expanding from three elements per side to six elements per side. 
This subdomain has proven very useful for flows having small areas of high 
velocity gradient where a very fine mesh is required relative to the overall 
flowfield mesh size, and windowing-in is not practical. Fig. A. 5 illustrates 
such a case where vortex generation at a wing tip is to be computed. Very 
high velocity gradients occur in the immediate area of the tip and decay 
rapidly. Utilizing the telescoping subdomain, the grid size is reduced from 
six elements per side in the tip region to two elements per side where an 
infinite boundary condition was applied and 64 elements were eliminated from 
the computational domain. 

The interpolation functions described above are used to distribute the 
gridpoint data over each subdomain. Combining the subdomains to form the 
complete discretization, however, poses the problem of having duplicate sets 
of data along the interconnecting boundaries. Removal of the duplication is 
accomplished through the construction of a subdomain connection table from 
the specified data and interrogating it from an algorithm which operates over 
the subdomain boundaries. Table A.l illustrates the connection table for the 
first four sub-domains of Fig. A.l. The storage requirements for the table 
are 2£(no. of subdomain sides). 
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TABLE A.l Subdomain Connection Table for Joukowski Airfoil. (First 4 Subdomains ) 


Subdomain 

Side 

Subdomain 

Side 

1 

1 

10 

3 

1 

2 

0 

0 

1 

3 

2 

1 

1 

4 

0 

0 

2 

1 

1 

3 

2 

2 

0 

0 

2 

3 

3 

1 

2 

4 

0 

0 

3 

1 

2 

3 

3 

2 

0 

0 

3 

3 

4 

1 

3 

4 

0 

0 

4 

1 

3 

3 

4 

2 

0 

0 

4 

3 

5 

1 

4 

4 

0 

o 


Note that the first two columns need not be stored since they are sequential. 

For the solution domain of Fig. A.l, therefore, only eighty words of storage 
are required to administrate boundary duplication elimination. The algorithm 
proceeds subdomain by subdomain, first generating a sequence of local 
(dummy) node numbers, then interrogating the connection table to determine 
if any sides are connected to any other subdomain. Upon finding one which 
is of lower number it retrieves the local boundary node numbers for that sub- 
domain and substitutes them for the coincident set into the proper locations. 
Simultaneously, the duplicate generated variable set at the boundary is elim- 
inated. The finite element connection table is subsequently formed for all 
generated elements in the subdomain and sequential global node numbering is 
substituted in a simple loop at the end. The connection table is also use- 
ful for locating external boundaries since they appear as zeros in the table. 
Note in the above scheme that only the sub-domain data was manipulated to 
sequence the generated data thus maintaining higher efficiency of the algorithm. 
Also, storage requirements are minimized by requiring operation on not more 
than two subdomains simultaneously, hence the data for a large problem could 
be generated interact! - vely on a minicomputer or small time sharing system 
utilizing disk or tape intermediate storage. An especially useful feature of 
such a method when combined with a graphics package and a video (CRT) terminal 
is ability to interactively modify and debug data decks. The above scheme is 
very efficient and generated data for 24 variables over 10 nodes/CPU sec on 
an IBM 360-65- 

The methodology described above requires a problem description made up 
of gridpoint coordinate and parametric data and also a list of gridpoint 
numbers on each subdomain. It has been found useful in aerodynamics to further 
automate data specification due to the overall similarity of the airfoil flow 
domains encountered. The flow domain of interest for single element airfoils 
is generally as shown in Fig. A.l with variations in shape of the boundary 
which is small compared to the overall flow domain size. The subdomain data, 
therefore, is program set except for specification of airfoil and wake 
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coordinates. These coordinates' are input in the form of a user supplied 
subroutine and interpolated using cubic splines to obtain coordinates at the 
subregion vertices The side nodes between the airfoil and in infinite boundary 
are automatically relocated as a function of the airfoil thickness to provide 
a consistant element size along the airfoil surface. Where the subroutine 
can be made sufficiently general to include coordinates as a function of maxi- 
mum thickness ratio, etc., such as for the Karmann-Trefftz series, a wide range 
of shapes may be analyzed by simply respecifying the key parameters such as 
thickness ratio, camber, trailing edge thickness, etc. Iterative viscous- 
inviscid solutions are achieved by automatically adding the boundary layer and 
wake thickness from the viscous solution to the basic airfoil shape and pro- 
ceeding through a new C p solution. Details on using the discretizer are given 
in the computer program^ users 1 manual (ref. 23). 


83 



REFERENCES 


1. Karemcheti , K. , Principles of Ideal -Fluid Aerodynamics , John Wiley, 

New York, 1966. 

2. Chapman, D.R., Mark, H., and Pirtle, M.W. , "Computers vs. Wind Tunnels," 

AIAA A/A. , 1975, pp. 22-30. 

3. Hess, J.L. and Smith, A.M.O., "Calculation of Potential Flow About 

Arbitrary Bodies," in Progress in Aero. Sciences , V. 8, PergamoTi 
Press, N.Y. , 1966. 

4. "Aerodynamic Analyses Requiring Advanced Computers," 

NASA SP -347, Parts I & II, 1975. 

5. "Vortex Lattice Utilization," NASA SP-405, 1976. 

6. Schlichting, H., Boundary- Layer Theory , McGraw-Hill, New York, 1968. 

7. Stevens, W.A., Goradia, S.H., and Braden, J.A., "Mathematical Model 

for Two-Dimensional Multi-Component Airfoils in Viscous Flow," 

NASA CR-1843, 1971. 

8. Zienkiewicz, O.C., The Finite Element Method in Engineering Science , 

McGraw Hill, London, 1971. 

9. Zienkiewicz, O.C. and Cheung, Y.K., "Finite Elements in the Solution 

of Field Problems," The Engineer, 1965, p. 507-510. 

10. deVries, G. and Norrie, D.A., "The Application of the Finite Element 

Technique to Potential Flow Problems," Trans. ASME, J. App. Mech., 
1971, p. 798-802. 

11. Meissner, Udo, "A Mixed Finite Element Model for Use in Potential Flow 

Problems." Int. J. Num. Mtd. Engr., V.6, p. 467-473, 1973. 

12. Isaacs, Lewis T. , "A Curved Cubic Triangular Finite Element for Potential 

Flow Problems," Int. J. Num. Mtd. Engr., V.7, p. 337-344, 1973. 

13. Vooren, J.V.D. and Laburjere, Th. E., "Finite Element Solution of the 

Incompressible Flow Over an Airfoil in a Nonuniform Stream," 
Proceedings Int. Conf. on Num. Mtd. in Fluid Dynamics, Southampton 
England, Sept. 1973. 

14. Hirsch, C. and Warzee, G., "A Finite Element Method for Flow Calculations 

in Turbo - Machines," Free University of Brussels Report V.U.B-Str-5, 
July 1974. 

15. Bratanow, T. and Ecer, A., "On the Applications of the Finite Element 

Method in Unsteady Aerodynamics," AIAA J., 12_, No. 4, p. 503-510, 

1974. 


84 


16. Ballhaus, W.F., Bailey, F.R., and Frick, J., "Improved Computational 

Treatment of Transonic Flow About Swept Wings," NASA CP-2001, 

Vol . 4, pp. 1311-1320, 1976, 

17. Cebeci , T. and Smith, A.M.O., Analysis of Turbulent Boundary Lavers , 

Academic Press, New York, 1974. 

18. Tennekes, H. and Lumley, J.L., A First Course in Turbulence , The MIT 

Press, Cambridge, USA, 1974. 

19. Launder, B.E., Reece, G.J. and Rodi , W., "Progress in the Development 

of a Reynolds-Stress Turbulence Closure," J. Flu. Mech., Vol. 68, 
Pt. 3, pp. 537-566, 1975. 

20. Hanjalic', K. and Launder, B.E., "A Reynolds Stress Model of Turbulence 

and its Application to Thin Shear Flows," J. Flu. Mech. Vol. 52, 
pp. 609-638, 1972. 

21. Launder, B.E. and Spalding, D.B., Lectures in Mathematical Models of 

Turbulence , Academic Press, London, 1972. 

22. Patankar, S.V. and Spalding, D.B., Heat and Mass Transfer in Boundary 

Layers , Second Ed., Intertext, London, 1970. 

23. Manhardt, P.D., Orzechowski, J.A., and Baker, A.J., "C0M0C II: Two- 

Dimensional Aerodynamics Sequence, Computer Program User's 
Guide," NASA CR-145165, 1977. 

24. Strang, G., Fix, G.J., An Analysis of the Finite Element Method , 

Prentice Hall, New Jersey, 1973. 

25. Shames, Irving H., Mechanics of Fluids , McGraw Hill, New York, 1962. 

26. Proceedings, AFOSR-IFP-Stanford Conference on Computation of Turbulent 

Boundary Layers - 1968, Vol. I., Eds. S.J. Kline, G. Sovran, M.V. 
Morkovan, D.J. Cockrell, Vol. II, Eds. D.A. Coles, E.A. Hirst, 
pub. by Thermosciences Div., Dept. Mech. Engr., Stanford University 
1969. 

27. Smith, A.M.O., and Cebeci, T., "Remarks on Methods for Predicting 

Viscous Drag," AGARD-CP-124, 1973. 

28. Preston, J.H. and Sweeting, N.E.,"The Experimental Determination of the 

Boundary Layer and Wake Characteristics of A Simple Joukowski 
Airfoil With Particular Reference to the Trailing Edge Region," 
Reports and Memoranda No. 1998, N.P.L., 1943 

29. Goradia, S.H., and Lilley, D.E., "Theoretical and Experimental Study 

of A New Method for Prediction of Profile Drag of Airfoil Sections, 
NASA CR-2539 , 1975. 


NASA-Langley, 1977 CR-2908 


85 



